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FINITE  ELEMENT  METHODS 

FOR  SINGULAR  TWO-POINT  BOUNDARY  VALUE  PROBLEMS 


Robert  Samuel  Schreiber 
Yale  University,  1977 


Ktntiii'n  etin 


Id)  a 

/l_J_ 


Until  quite  recently,  few  effective  numerical  solution  techniques 
were  known  for  solving  two-point  boundary  value  problems  for  the 
equation 


-;^,(p(x)^)  + qu  - f,  0<x<l,  p(0)  - 0. 


'dx' 


In  this  dissertation  we  analyze  several  new  finite  element  methods  for 
approximating  the  solution  of  this  problem,  and  present)  new  analyses  for 
some  known  methods.  , 


Two  classes  of  singular  problem  may  be  distinguished,  depending  on 

whether  or  not  p(x)  ^ l\o,1).  In  the  first  of  these,  p(x)  behaves 

o 1-0 

like  X , 0£O<l,  near  0.  The  solution  u(x)  has  a singularity  (like  x ) 

at  0,  so  that  its  derivatives  there  are  infinite.  This  frustrates  all 

the  usual  theory.  The  finite  element  techniques  proposed  so  far  require 

chat  the  basis  Include  functions  which  mimic  the  behavior  of  the 

solution. 


We  investigate  the  idea  of  approximating  Che  solution  with 
piecewise  polynomials  on  a nonuniform  mesh  adapted  to  the  singularity. 
Given  knowledge  of  the  singularity  (i.e.,  o),  it  is  possible  to 


construct  a sequence  of  graded  meshes  such  that  the  rate  of 


2 

convergence  in  the  L -norm  is  the  best  possible.  We  prove  upper  and 
lower  bounds  on  the  extent  to  which  the  mesh  must  be  graded. 

Tlie  solution  can  also  be  approximated  by  a function  of  the  form 

-0 

X s(x),  with  s(x)  a piecewise  polynomial.  We  obtain  error  bounds  and 
numerical  results  for  these  "weighted  splines"  which  indicate  that  they 
are  the  best  for  practical  computation.  For  a third  subspace  (due  to 
Crouzeix  and  Thomas),  we  Improve  known  error  bounds  by  using  a mildly 
graded  mesh. 

Problems  of  the  second  kind,  where  p(x)  behaves  like  x near  0, 
with  k a positive  integer,  arise  from  spherically  symmetric  elliptic 
boundary  value  problems  in  n **  k+1  dimensions.  The  solutions  are 
smooth,  indicating  that  piecewise  polynomials  on  uniform  or  nearly 
uniform  mesh  will  be  effective  approximators. 

We  improve  previous  results  by  removing  any  restriction  on  n,  by 
adding  the  requirement  that  the  approximating  spline  functions  be  smooth 
at  the  center  of  the  n-dimensional  domain,  and  by  obtaining  error  bounds 
(of  the  best  possible  order)  in  the  natural  norms  for  the  problem. 
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CHAPTER  1 
INTRUDUCTIUN 


i.  1 ^ ilx ample 

Consider  the  two-point  boundary  value  problem 

(1.1)  — “ "j*  0 < X < 1, 

with  boundary  conditions 

(1.2)  u(0)  - u(l)  - 0. 

1/2  3/2 

The  unique  solution,  u(x)  - x - x , is  (uniformly)  continuous  on 
[0,1]  and  continuously  differentiable  (indeed  analytic)  in  (0,1). 

However,  its  derivative  grows  arbitrarily  large  as  x approaches  0; 

2 2 
consequently,  while  u € L (0,1),  none  of  its  derivatives  is  in  L (0,1). 

Therefore,  analyses  of  numerical  solution  techniques  for  two-point 

boundary  value  problems  which  require  that  certain  derivatives  of  u be 

bounded  don't  apply  to  this  problem. 

Equation  (1.1)  falls  outside  the  realm  of  problems  treated  in  most 
discussions  of  numerical  methods  for  boundary  value  problems  (both  for 


- 2 - 


ordinary  and  elliptic  partial  differential  equations)  because  the 
coefficient  /x  vanishes  at  a point  of  the  boundary  of  the  domain.  This 
"singular  point"  Is  precisely  where  the  derivatives  of  the  solution  blow 
up,  making  It  Impossible  to  apply  well-known  numerical  techniques  to  the 
problem.  The  goal  of  this  research  Is  to  extend  numerical  methods, 
specifically  the  finite  element  method,  to  apply  to  singular  two-point 
boundary  value  problems.  Ue  hope  to  develop  new  methods  for 
approximating  the  solution  and  tc  justify  those  methods  with  rigorous 
error  bounds. 

1 . 2 Summary  of  Results  Obtained 


We  shall  treat  equations  of  the  form 


(1.3) 


- + q(x)u  - f{x),  0 < X < 1, 


where  p(x)  > 0 for  x e (0,1],  and  p(0)  - 0.  (In  the  nonsingular  case, 

it  Is  assumed  that  p(x)  > 0 for  all  x e [0,1].)  Two  classes  of  singular 

problem  may  be  distinguished,  depending  on  whether  or  not  the  Integral 
1 

/ is  convergent,  i.e.,  finite. 


Part  II  of  the  dissertation  Is  concerned  with  the  first  of  these 
two  types  of  singular  two-point  boundary  value  problem.  Such  problems 
arise  In  potential  theory.  A function  u which  satisfies  the  equation 


L„lu] 


3^u  . a^u  . a au 

71  71 

3x^  ay^ 


is  said  to  be  a generalized  axially  symmetric  potential  in  n » o+2 
dimensions;  o need  not  be  integral  (see  [Wl]).  Parter,  in  his  early 
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treatment  of  numerical  methods  for  generalized  axially  symmetric 
potentials  In  a rectangle,  arrived  at  the  equation 

v(0)-l,  v(l)-0,  |ol<l, 

by  separation  of  variables  [PI].  Jamet  [Jl]  later  considered  the  more 
general  problem 

(1.4)  - D(x'^p(x)Ou)  + q(x)u  “ f(x),  0<x<i,  0£o<l, 

with  the  boundary  conditions  (1.2),  where  p Is  a smooth  function 
strictly  positive  in  the  Interval  [0,1].  For  a more  complete  account  of 

i 

earlier  work  on  numerical  methods  for  this  problem,  see  Section  4.1.  j 

The  difficulty  with  these  problems  Is  that  their  solutions  are  not 
smooth  In  the  usual  sense  of  having  two  more  derivatives  than  f;  In  ' 

fact,  their  derivatives  are  unbounded  at  the  origin.  Nevertheless  we 
are  able  to  define  and  analyze  two  new  numerical  approximation  schemes 
of  hlgh-order  accuracy. 

As  approximators  of  the  solutions  of  (1.3),  we  shall  use  functions 
which  are  splines  (piecewise  polynomials)  or  are  in  some  way  related  to 
such  functions.  A function  is  a spline  on,  say  [0,1],  If  It  is  a 
polynomial  of  some  fixed  degree  on  each  Interval  of  a given  partition  of 
[0,1].  In  addition,  various  continuity  requirements  may  be  imposed; 
thus  we  may  speak  of  C^-cubics  — functions  which  are  cubic  polynomials 
in  each  interval  and,  together  with  their  first  derivatives,  are 
continuous  at  the  "knots"  of  the  partition. 


J 


The  essential  point  about  the  solution  u of  (1.4)  - (1.2)  Is  that 
u(x)  behaves  like  ° near  0.  (Thus,  for  example,  v(x)  = x°u(x)  Is 
smooth  and  vanishes  at  0 (Theorem  3.1)).  This  leads  Immediately  to  the 
following  Idea.  It  has  been  shown  that  piecewise  polynomials  can 
approximate  nonsraooth  functions  of  the  form  x**,  o not  an  Integer, 
essentially  as  well  as  smooth  functions,  simply  by  using  a nonuniform 
mesh.  Thus,  It  should  be  possible  to  approximate  the  solution  of 
(1.4)  - (1.2)  with  piecewise  polynomials  on  a nonuniform  mesh  adapted  to 
the  singularity.  We  show  how,  given  knowledge  of  the  singularity  (l.e., 
o) , to  construct  a sequence  of  meshes  such  that  the  rate  of  convergence 
Is  the  best  possible. 

As  an  aside,  we  give  a problem  and  a subspace  (one  Involving  a 
partition  which  Is  not  sufficiently  skewed  towards  0)  for  which  the  RRG 

approximation  does  not  converge  at  the  best  possible  rate  In  the 

2 2 
L -norm;  It  Is  worse  by  an  order  of  magnitude  than  the  L -projection  of 

the  solution.  This  appears  to  be  the  first  such  example  reported.  It 

has  been  shown  (Elsenstat,  Schrelber,  and  Schultz  [El])  that  this  cannot 

happen  for  a wide  class  of  nonslngular  problems. 


1 


We  next  considet*  subspaces  of  "weighted  splines,"  functions  of  the 

. 0 _ 0 

form  X s(x)  with  s(x)  a spline.  Since  we  Know  that  u(x)  “ x v(x) 

where  v Is  smooth.  It  seems  likely  that  weighted  splines  will  be  good 
approximators  of  u.  We  obtain  error  bounds  which  show  that  this  is 
Indeed  the  case.  These  spaces  appear  to  be  the  best  for  practical 


computation. 


I 

i 


f 
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Work  on  this  problem  using  the  finite  element  method  was  begun  by 
Ciarlet,  Nattercr,  and  Varga  [Cl],  who  used  subspaces  of  functions  which 
are  piecewise  elements  of  the  null  space  of  the  differential  operator. 
These  "L-spllnes"  were  generalized  by  Crouzeix  and  Thomas  [C4] , who 
investigated  functions  which  are  mapped  into  polynomials  (not  just  the 
zero  polynomial)  by  the  operator.  Applying  these  spaces  to  the  problem 

(1.4) ,  Crouzeix  and  Thomas  obtained  energy-norm  and  L -norm  error  bounds 
which  are  high-order,  but  not  quite  optimal.  We  present  their  theory 
(with  some  simplifications)  in  Chapter  8.  We  then  improve  their  results 
by  showing  that  if  a mildly  graded  partition  is  used,  then  the  error  in 
their  procedure  is  of  optimal  order. 

Part  111  of  the  dissertation  is  concerned  with  singular  problems 

k 

(1.3)  in  which  p(x)  = x where  k is  a positive  Integer.  The  boundary 
conditions  are  not  those  of  (1.2),  but  rather, 

(1.5)  u(l)  = 0,  Du(0)  = 0 (alternatively,  u(0)  finite). 

These  problems  arise  from  multidimensional  elliptic  boundary  value 
problems  (for  example,  the  Dirichlet  problem  for  Poisson's  equation 
-AU  = F)  which  possess  spherical  symmetry.  The  solution  is  a function 
only  of  distance  r from  the  origin  in  IR^  and  can  be  obtained  by  solving 
a singular  two-point  boundary  value  problem  of  this  type.  Unlike 
problems  of  the  first  type,  the  solution  is  smooth.  Nevertheless,  the 
usual  theory  for  two-point  boundary  value  problems  breaks  down  when 
applied  here.  A discussion  of  earlier  work  on  numerical  methods  for 


( 


I 

I 

\ 

i 

i 

I 

! 


i 

\ 


this  problem  is  given  in  Section  9.1. 
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We  propose  to  approximate  the  solution  of  such  an  equation  by 
spherically  symmetric  functions  which  are  splines  in  the  variable  r 
(distance  from  the  center  of  the  domain).  This  procedure  has  been 
analyzed  previously  ( [D6]  , (J3])  for  the  special  cases  n ■ 2 and  3,  with 
error  bounds  obtained  in  the  usual  Sobolev  norms  on  the  Interval  [0,1). 
We  improve  these  results  in  three  directions:  by  removing  any 
restriction  on  n,  by  adding  the  requirement  that  the  approximating 
spline  functions  be  smooth  at  the  center  of  the  n-dlmenslonal  domain, 
and  by  obtaining  error  bounds  (of  the  best  possible  order)  in  the 
"natural"  norms  for  the  problem  — the  Sobolev  norms  on  the  original 
domain  in  IR”  Instead  of  the  interval  to  which  the  problem  has  been 
reduced. 


1 . 3 Outline  of  the  Dissertation 

There  are  three  main  parts  to  the  dissertation,  the  first  of  which 
is  devoted  to  mathematical  preliminaries.  The  numerical  methods  we 
employ  are  all  cases  of  a general  class  of  approximation  techniques  for 
solving  boundary  value  problems:  the  Raylelgh-Rltz-Galerkin  (RRG) 
method.  Chapter  2 is  a discussion  of  the  Rayleigh-Ritz-Galerkin 
procedure  as  applied  to  singular  two-point  boundary  value  problems. 
Chapter  3 introduces  the  spaces  of  spline  functions  to  be  used  as 
approximate  solutions. 

In  Part  II,  we  consider  the  problem  (1.4)  - (1.2).  In  Chapter  4, 


we  introduce  the  variational  form  of  the  problem  in  an  appropr  '.ate  class 
of  function  spaces.  One  of  the  underlying  ideas  of  Part  II  is  the  use 


7 


of  splines  on  a class  of  graded  nonunlfom  partitions  of  [0,1].  These 
"B“graded"  meshes  have  previously  been  shown  to  be  useful  in 
approximating  functions  with  a singularity.  We  define  the  B-graded 
meshes  and  summarize  several  of  their  properties  in  Section  4.3. 
Chapter  5 is  a compendium  of  results  on  the  properties  of  the  solution 
which  will  be  used  throughout  Part  11.  In  Chapter  6,  we  consider  the 
approximation  of  singular  functions  (elements  of  tlie  spaces  introduced 
in  Chapter  4)  by  piecewise  polynomials  on  a B-graded  mesh,  and  in 
Chapter  7,  we  consider  approximation  by  weighted  splines.  Numerical 
results  are  included  in  both  chapters.  In  Chapter  8,  we  deal  with  the 
generalized  L-splines  of  Crouzelx  and  Thomas. 

Part  III  of  the  dissertation  is  concerned  with  singular  problems 
(1.3)  - (1.5)  in  which  p(x)  behaves  like  x near  0,  where  k is  a 
positive  integer.  In  Chapter  9,  we  obtain  the  variational  form  of  the 
problem,  and  in  Chapter  10,  we  define  and  analyze  several  spline 
approximation  schemes.  Numerical  results  are  included. 

I 

I 


I 
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CHAPTER  2 


THE  RAYLEIGH-RITZ-GALERKIN  PROCEDURE 
FOR  SINGULAR  TWO-POINT  BOUNDARY  VALUE  PROBLEMS 


2. 1 Notation 

In  this  section  we  Introduce  notation,  definitions,  and  results 
which  will  be  used  throughout  this  dissertation. 

Let  1 = [0,1].  For  each  Integer  t > 0,  let  D*'f(x)  = ^(x).  For  S 

dx^ 

a finite  set,  let  |S|  denote  the  number  of  elements  of  S. 

All  the  functions  we  deal  with  will  be  real-valued.  The  support  of 
a function  f,  supp(f),  is  the  closure  of  the  set  of  points  x such  that 
I f(x)  0.  If  f is  a bounded  function  on  (a,b),  we  write 


l|£  II 


L"(a,b) 


sup  I f (x) I . 
a<x<b 


For  m a non-negative  integer,  let  H™(a,b)  (respectively,  H®(a,b)) 
b3  the  closure  of  the  C"  functions  (respectively,  the  C“  functions  with 
compact  support  in  (a,b))  with  respect  to  the  norm 


8 


H"'(a,b) 


f 
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I ^ (D-)f(x))2  dx 
j“0  a 


These  are  known  as  Sobolev  spaces.  Note  that  H™(a,b)  may  be  Identified 
with  the  space  of  C*”  ^ functions  with  absolutely  continuous 
(tn-l  )*'^-derivatlve  and  m*'^'-der ivativc  in  L^(a,b),  while  HQ(a,b)  may  be 
Identified  with  the  subspace  of  H^Cajh)  consisting  of  functions  which 
vanish,  along  with  their  first  m-1  derivatives,  at  a and  b.  Thus, 
H^(a,b)  = L^(a,b).  H'"(a,b)  is  a Hilbert  space  with  inner  product 


(f,g) 


H"'{a,b) 


mb 

Z / f D-*  g dx . 


j=*0  a 

For  a complete  discussion  of  Sobolev  spaces,  see  Adams  (All. 

We  also  define  VT ’ (a,b)  (respectively,  W™’  ) as  the  closure  of 
C°°(a,b)  (respectively,  {f  t C”(a,b)  | f has  compact  support  in  (a,b)  }) 
with  respect  to  the  norm 


m 

r 


(a,b)  j=0 


D-^f 


L"(a,b)’ 


We  denote  by  h’"  (respectively,  Hq,  L°°,  W™’  ) the  space  H™(0,1) 

00 

(respectively,  Hq(0,1),  L°°(0,1),  w"*’  (0,1)).  Moreover,  we  denote  the 


norm  and  inner  product  of  by 


and  ( . , . ) , and  let 
m ra 


L“’(0,1)' 


Theorem  2. 1 (Rayleigh-Ri tz  Inequality)  [Bl]:  If  f e HQ(a,b),  then 

9^0  9^9 

/ f dx  < (b  - a)  -C  (Df)‘‘  dx. 
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Theorem  2.2  (Markov's  Inequality)  (Ml):  For  all  polynomials  p of 

n 


degree  n and  all  real  a,  b with  a < b. 


"^n‘'L"(a,b)  1 " (b  - a)'  II  P Jl  L-’Ca.b) ' 


In  Part  111,  we  will  be  concerned  with  functions  defined  on  a 
bounded,  open  subset  B of  IR".  We  define  the  space  H®{B)  by  taking  the 
closure  of  the  C°°-f unctions  with  respect  to  the  norm 


( 


h“(b) 


I f (D“f)^  dx 
|a|<m  B 


where 


|a| 


V.  V 

1 z n 


Sobolev's  lemma  [F3,  p.  283]  Is  concerned  with  bounds  on  the  value 

2 

of  a function  at  an  arbitrary  point  In  terms  of  L -norms  of  its 
derivatives.  While  the  following  theorem  can  be  proved  for  a more 
general  class  of  domains,  we  do  not  require  anything  more  than  a ball  in 
IR”  (an  Interval  In  IR).  Let  L^J  denote  the  largest  integer  not 
exceeding  x. 

Theorem  2.3;  Let  B be  a ball  In  IR".  If  u e H®(B) , where 
m = LtJ  + 1,  then  u can  be  identified  with  a uniformly  continuous 
function  u(x)  In  B such  that 


^ lli^fR^  1 C Hu  11 

L (B)  - n ^m^gj 


where  is  a constant  which  depends  only  on  n. 


w 
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In  the  one-dlmenslonal  case,  Sobolev's  lemma  shows  Chat 
H^(a,b)  c L"(a,b)  and 

^ L"(a,b)  - ^1  II  ^ 'I  H‘  (a.b)  ^ ^ ” (a,b) . 


If,  In  addition,  f(a)  = 0,  we  have  a stronger  bound. 


Lemma  2. 1:  If  f c H^(a,b)  and  f(a)  = 0,  then 


(2.1)  II  ^ II  L^fa.b)  - (b  - a)  ^ H H L*(a,b)  * 


Proof:  Since  f is  absolutely  continuous. 


f(x) 


; Df(t)  dt 
a 


Using  the  Cauchy-Schwarz  Inequality, 


lf(x)  I < (x  - a) 


1/2 


f (Df(t))2  dt 


< (b  - II  Df 


L"(a,b)- 


We  will  sometimes  use  the  letter  C to  denote  a generic  constant, 
not  the  same  at  each  occurrence.  Within  each  chapter  we  define  "local" 
constants  Cq,  Cp  ...,  to  be  used  and  referred  to  only  within  that 


chapter 
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2 . 2 The  Raylelgh-Rltz-Galerkln  Procedure 
We  consider  Che  equation 

(2.2)  Lu  - f, 

where  L Is  a densely  defined  linear  operator  on  a (real)  Hilbert  space  H 

and  f e H.  The  development  follows  Mlkhlln  [M2],  to  which  the  reader  Is 

referred  for  further  details  and  proofs.  In  the  sequel,  we  shall  be 

concerned  with  Hilbert  spaces  of  functions  defined  on  a real  Interval 

2 

(a,b).  In  Part  11,  L (a,b)  plays  the  role  of  H,  while  In  Part  III, 
where  we  study  elliptic  problems  In  a ball  In  Ir",  H Is  the  space  of 
functions  square  Integrable  on  (a,b)  with  respect  to  the  weight  function 
x**  L Is  a symmetric  differential  operator  with  principal  part 

-D(p(x)Du)  and  p(0)  =•  0.  | 

i 

We  assume  that  L Is  defined  on  a dense  subspace  M of  H and  that  L j 

Is  symmetric  and  positive  definite:  for  every  u,v  £ M, 

(2.3)  a(u,v)  = (Lu,v)  « (u,Lv)  | 

and  I 

(2.4)  ||v||2  H a(v,v)  > i^rllvllj,  ] 

E H 1 

where  6 Is  a constant  Independent  of  v.  It  Is  well-known  that  the  \ 

closure  S of  M with  respect  to  the  norm  II  . II  ^ is  a Hilbert  space  (and  a 

I 

subspace  of  H)  with  inner  product  a(u,v)  and  corresponding  norm  ||  • II  j 

4 

i 

Theorem  2.4  [S4]:  For  every  f e H,  there  exists  a unique  element  u e S, 

called  the  generalized  solution  of  (2.2),  satisfying  i 

J 

i 


t 
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(2.5)  a(u,v)  - (f.v)jj  for  all  v e S. 

Moreover,  u e S satisfies  (2.5)  if  and  only  if  it  minimizes  the 
quadratic  functional 

F[u]  = a(u,u)  - 2(f,u)^,  u € S. 

The  solution  u is  bounded  in  terms  of  the  data  f.  For,  by  (2.5) 
and  the  Cauchy-Schwarz  inequality, 

II  u II  ^ - a(u,u)  = 

< II  f II  H II  u II  H 

< II  f II  H 6 II  u II  £, 

whence 

(2.6)  Hu  11^  < i\\i\\^, 

and  by  (2.4), 

(2.7)  llulljj  < S^llflljj. 

Let  be  a finite-dimensional  subspace  of  S.  The  Galerkin 
approximation  to  u in  Is  the  unique  element  u of  satisfying 

(2.8)  a(u.Vj^)  = (f,v^)jj  for  all  e Sj^. 

The  Raylelgh-Ritz  approximation  is  the  unique  element  of  which 
minimizes  the  functional  F[u]  over  as  in  the  case  of  the 

generalized  solution,  these  definitions  are  equivalent;  we  call  this 
element  the  Rayleigh-Ritz-Galerkin  (RRG)  approximation  to  u in  S . 


u - 


Theorem  2 . 5;  The  Raylelgh-Rlcz-Calerkln  approximation  u is  the  best 
possible  approximation  to  u in  the  subspace  S with  respect  to  the  norm 

II  * II  £,  i ae  • , 

(2.9)  II  u - u II  ^ . inf  II  u - (( 

V eS 
n n 

Proof ; Subtracting  (2.8)  from  (2.5)  shows  that  the  error  u - u is 
orthogonal  to  S^,  i.e., 

a(u  - u,v  ) = 0 for  all  v„  c s„. 

n n n 

The  bound  (2.9)  then  follows  by  a standard  argument  (see  [S31). 

a 


Suppose  we  want  to  determine  the  coefficients  of  the  RRG 

approximation  with  respect  to  a basis  (B.,  B } for  S . Bv 

(2.8), 


a(u,B  ) - (f,B.)jj,  1 < j < n. 


Thus,  if 


a - ^ ^iB,, 


i-i 


then  ^ is  the  solution  of  the  linear  system  of  equations 


(2.10)  Ai  - f.  A - [a(B^,Bj)],  f-  [(f,B^)jj]. 


I 


J 
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The  equations  (2.10)  have  a unique  solution,  since  the  matrix  A Is 
symmetric  and  positive  definite.  In  fact.  Its  symmetry  Is  apparent,  and 
If  n - (n^,  02,  ....  n^)'^  c m",  then 


'j'  n n 

n Aji  - a(  £ n,B4,  ^ n.Bi) 

1-1  1-1 


i OiBi  II  I > 0 


with  equality  If  and  only  If  - 0.  We  will  later  Introduce  subspaces 
with  computationally  attractive  local  bases  (the  functions  are  zero 
over  most  of  (a,b)).  The  resulting  linear  systems  will  be  sparse  and, 
for  moderate  sized  systems,  we 11- conditioned. 


When  H Is  a function  space  on  (a,b)  and  L Is  a differential 
operator  of  the  form  (1.3),  the  coefficients  p and  q of  L enter  Into  the 
energy  norm,  which  makes  It  cumbersome  to  use  In  approximation- 
theoretic  arguments.  Therefore,  we  will  work  with  an  equivalent  norm 


g;  we  assume  that 


(2.11)  ^l|v||  < ||v||  £ ^l|v||g  forallveS, 


where  X and  A positive  constants  Independent  of  v.  (In  particular. 


/ w(x)(Dv)^  dx. 


where  in  Part  II,  w(x)  = x*’,  and  in  Part  III,  w(x)  = x"”^).  The  a 
priori  bound 


(2.12) 


then  follows  from  (2.6);  by  (2.9)  we  obtain  the  error  bound 
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(2.13) 


II  u - u II  s < y 11  U - U 11  g 

< { ln£  Hu-  11  £ 

< j int  Hu-  V 11 

n n 


2.3  Error  Bounds  for  the  RRG  Approximation 

In  this  section,  we  think  of  the  subspaces  as  being  piecewise 
polynomial  spaces  with  respect  to  a mesh  of  N sub  Intervals.  For  such 
spaces,  the  dimension  n Is  a linear  polynomial  In  N. 

The  RRG  approximation  u Is  the  best  possible  energy-norm 
approximation  to  u in  the  subspace  (equation  (2.9)).  Therefore,  the 
first  step  in  bounding  u - u Is  to  find  upper  bounds  for 
inf  II  u - V II  r*  This  is  generally  done  by  defining  an  approximation 

n bt 

V es 
n n 

mapping  from  S into  S (Interpolation  is  a typical  example)  for  which 

n 

error  bounds  can  be  obtained  a priori.  Any  such  error  bound  Implies  a 
corresponding  bound  for  ||  u - u II  j.. 

By  (2.4),  any  bound  on  ||  u - u II  automatically  induces  a bound  on 

E 

II  u - u II  „.  However,  the  resulting  bound  is  generally  not  sharp,  in  the 
sense  that  the  dependence  on  N (the  order  of  the  approximation)  is  not 
the  best  possible.  When  optimal-order  error  bounds  in  the  H-norm  are 
desired,  more  machinery  must  be  developed.  The  following  argument 


(Nltsche's  trick)  was  first  used  by  Nitsche  [N2]. 
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We  make  the  following  "approximation  hypothesis"  concerning  the 
ability  of  elements  of  to  approximate  solutions  of  the  problem  (2.2). 

(Al)  There  exists  an  integer  k > 2,  subspaces 

\ X2  ^ S,  and  a constant  Aj  such  that,  if  v ^ X^, 

then 

inf  llv  - V II  < A.  IIvII  . 

V es  n & 1 

n n 

In  Part  II,  we  take  II  u ||  ^ =•  ||D(x°Du)  II  j_2,  while  in  Part  III, 

II  “ II  X “ 11^^  just  the  usual  L^-norm  of  D^u  over 

the  unit  ball  in  IR”. 

We  also  need  a "regularity  hypothesis",  giving  bounds  on  the 

X^-norm  of  a generalized  solution  of  (2.2)  in  terms  of  the  H-norra  of  the 
right-hand  side. 

(A2)  There  exists  a positive  constant  T such  that,  if  f £ H,  then 
the  generalized  solution  u £ X^  and 

Hu  llx^  - II  H- 

In  Part  II,  regularity  is  provided  by  Theorem  5.2;  in  Part  III,  by 


Lemma  9.2. 


r 
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Theorem  2 ■ 6 (NiCsche):  If  satisfies  the  approximation  hypothesis 
(Al),  the  problem  satisfies  the  regularity  hypothesis  (A2),  and  u e 
then 

(2.14)  Hu  - ullg  < Y Aj  IIu  IIx 

and 

(2.15)  llu  - u||  < (AA,)^r  n"^  IIu  ||„  . 

H 1 Xj^ 

Proof ; (2.14)  follows  Immediately  from  the  approximation  hypothesis  and 

the  S-norm  quasi-optimality  (2.13)  of  the  RRG  approximation. 

Let  't'  = u - u and  let  * be  the  generalized  solution  of  (2.2)  with 
right-hand  side  'f.  By  the  definition  (2.5)  of  the  generalized  solution, 

II  u - u II  y = ('*',u  - u)jj  = a('*’,u  - u). 

Let  ♦ e S be  a best  E-norm  approximation  to  Since  u - u is 
n n 

orthogonal  (in  the  a(.,.)  inner  product)  to  S , 

n 

a(^  ,u  - u)  = 0. 


Adding  the  two  previous  equations  and  using  the  Cauchy-Schwarz 
inequality. 
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By  the  approximation  hypothesis  and  the  optimality  of  the  RRC 
approximation  in  the  E-nonn  (2.9), 


u - u 


inf  II  u - V 


V eS 

n n 


n b 


< A inf  II  u - V 


V eS 
n n 


n " S 


< AA.  Hu  IIj^  . 


‘1 


Using  these  two  inequalities  to  bound  the  right-hand  side  of  (2.1b),  we 
obtain  (2. 15) . 


D 


CHAPTER  3 


ON  PIECEWISE  POLYNOMIALS 


[ 


3. 1 Introduction 

In  this  chapter,  we  introduce  a broad  class  of  piecewise  polynomial 
(spline)  subspaces,  including  those  used  in  practical  computation. 

After  defining  the  well-known  "B-spline"  basis  functions  and  stating 
several  of  their  properties,  we  discuss  approximation  by  splines, 
specifically,  the  quasiintef polant . The  main  result  is  a local  error 
bound  for  the  quasilnterpolant,  due  to  de  Boor  and  Fix  [D3).  We  also 
show  that  the  quasilnterpolant  can  be  made  to  interpolate  at  the  end 
points,  a fact  we  will  later  find  useful  in  showing  the  existence  of 
spline  approximations  satisfying  various  boundary  conditions. 


3.2  Piecewise  Polynomial  Approximation 

Let  A:a  = x^  < < . . . < Xj^  = b be  a partition  of  (a,b),  and 

define 


1 


1 


^^1-1 ’^1^ ’ 


1 < i < N, 


- 20 


1 
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■ ^1  - ^i-1  • 1 1 i 1 N, 

h = max  h 
l<i<N 

and  Che  local  mesh  ratio 

hi 

M(A)  = max  . 

Let  k be  a positive  integer  and  let  z - (z^,  z^,  the 

incidence  vector  associated  with  A,  have  positive  integer  components, 
earh  less  than  or  equal  to  k-1 ; i.e.,  1 _<  < k-1 , 1 < i < N-1. 


Definition  3.1:  The  real-valued  function  s(x)  is  a spline  of  order  k 
for  A and  ^ if  s(x)  coincides  with  a polynomial  of  degree  < k (i.e., 
decree  < k)  on  each  open  subinterval  of  A and  has  k-l-z^  continuous 
derivatives  at  each  x^,  i.e., 

D^sCx^-)  = D-^s(x^+),  0 _<  j < k-l-z^,  1 1 i < N-1. 

G 

The  class  of  all  splines  of  order  k for  A and  z is  denoted  by  S (A,^). 
k N-1 

The  dimension  of  S (A,^)  is  d = k + t z.. 

1=1  ^ 

T 

Let  ^ = (1,1,  ...,  1)  be  the  (N-1 )-dimensional  "one-vector".  We 
k V 

denote  by  S (A)  the  space  S (A,^)  and  call  this  the  space  of  smooth 
splines  of  order  k for  A.  If  z is  any  incidence  vector,  then 
S^(A)  c s'^lA,^).  We  denote  by  Sq(A^^)  (respectively,  Sq(A))  the 
subspace  of  functions  s e S^(A,^)  (respectively,  S^(a))  such  that 


s(a)  = s(b)  = 0. 


r 
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The  utility  of  splines  for  practical  computation  stems  largely  from 

the  existence  of  a we 11- conditioned  local  basis,  the  B-spllnes.  In 

order  to  define  the  B-splines,  we  need  one  additional  concept.  A vector 
d*fk 

^ ® k-extended  partition  of  (a,b)  provided 


a) 


S - 4 


(3.1) 


b) 


•^d+l  “ •^d+2 


“ '^d+k  “ 


< tj+i,  1 < j < d+k-l, 

“>  'j  ' tj.k.  I i j i d- 


Given  A and  let  _t(A,£)  be  the  k-extended  partition  such  that 
1)  tj  « I 0 < i < N}. 

il)  the  multiplicity  of  x^  in  { tj } is  z^. 

Let  f(xQ,  xj,  ....  xi^)  denote  the  k*"*^  divided  difference  of  the 
function  f on  the  points  ( Xq  , Xj , ...,  xj^}.  The  set  of  normalized 
B-splines  on  ^ is  defined  by 


“ (‘^j+k  - '^j>8k(Lj tj+k’  1 < j < d. 


where 


g.(s;  x)  » (s 


-x,r' 

[0  s < X 


We  note  the  following  properties  of  N , (x)  [D2]: 

J *k 
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a)  supp(Nj^^)  - l<j<d; 

I'^il  - l(j  I supp(Nj^j^)  n ^ (^}|  . k; 


(3*2)  c)  0 ^ N u(x)  £1,  a £ X < b; 


® J — 2+i. 


= 0 for  all  J £ d-1-1; 

e)  {Nj  I 1 £ j £ d}  is  a basis  for  S^(A,^), 


We  shall  require  a bound  on  the  derivatives  of  N 


Lemma  3. 1 : For  all  r £ 0 and  all  1 £ j £ d, 

0.3)  ID^N^  ^(x)  ( £ all  x e , 

where  B = n (k  - 1)^. 

" ^=l 

Proof : By  (3.2)(c),  (3.3)  holds  for  r = 0.  (We  adopt  the  convention 
0 

that  n I = 1.)  To  prove  the  general  result  (by  Induction  on  r) , we  use 
1=1 

Markov's  inequality  (Theorem  2.3).  Since  d’^n  , is  a polynomial  of 

J 

degree  k-r-l  in  each  of  the  intervals 

£ (k-r-l  )B^h]^^’^^^ 

- R ,-(>^1)  X I 

®r+l''l  • " ^i* 


G 
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Birkhoff  [B2]  devised  a projection  mapping  P from  spaces  of  smooth 
functions  into  S (a»£)  with  the  property  that,  in  the  interval  1^,  Pf 
depends  on  f only  in  some  small  neighborhood  of  and  the  rate  of 
convergence  of  Pf  to  f is  optimal.  His  work  was  generalized  to  n 
dimensions  (and  simplified  in  one  dimension)  by  de  Boor  and  Fix  [D3] , 
who  named  Pf  the  quasiinterpolant  of  f. 

Definition  3. 2 (de  Boor  and  Fix  [D3]):  For  each  Integer  j,  1 ^ j £ d, 
let  Tj  be  a point  in  the  support  of  Nj  l.e.,  tj  £ Tj  £ ^j+k'  ^ 

have  i-l  continuous  derivatives,  1 £ i,  £ k.  Then  the  quasiinterpolant 

(3.4)  F^£  , Xj£ 

where  the  linear  functionals  x,  are  given  by 


(3.5) 

(3.6) 
and 


X,f  = X.  f 
j J 


“3,£  » <-■> 


r<l 

k-l-r  , - 

k-l-r  ° *1  ^ ^ 


(k-D! 


(3.7)  it'j(x)  = - x)  . . . - x) 


) 


i 


I 

I 

I 

I 


w 
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Lemma  3.2  (de  Boor  and  Fix  [D3]):  If  p Is  a polynomial  of  degree  < I, 
then 


(3.8) 


F p = p. 


Let  be  the  smallest  Interval  containing  both  1^  and 
I j e v^}.  Clearly,  9^  ^*i-k'  *i+k-l^  Figure  3.1). 


In  I^,  F^f  depends  only  on  the  values  of  f in  0^;  thus  we  have  local 
error  bounds. 

I 

Theorem  3.  1 (de  Boor  and  Fix  (D3]):  Let  f £ H (0^).  There  exists  a 
positive  constant  K = K(k,  j ,M(  A) ) , such  that 

(3.9)  llDJ(f  - F^f)  ^ < K |0j^'J  IlD^f  IIl^O  ). 

where  |0^|  = length  of  9^. 

An  important  fact  about  the  quasiinterpolant  is  that  it  can  be  made 


to  interpolate  f,  and  its  derivatives  of  orders  1 through  1-1 , at  a 


J 
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and  b.  This  will  be  useful  in  proving  the  existence  of  spline  approxi- 
mations satisfying  various  boundary  conditions. 

Lemma  3. 3:  Let  the  first  (respectively,  last)  m £ quasiinterpolation 
points  be  placed  at  a (respectively,  b).  Then  the  quasiinterpolant 

F f interpolates  f and  its  first  m-1  derivatives  at  that  point. 

Proof : We  prove  the  result  for  interpolation  at  a;  the  result  for  b 

follows  by  symmetry.  According  to  (3.2)(d),  D^N . , (a)  =■  0 for  j > i+1, 

J 

so  that 


(3.13) 


. i+1 

D^F  f(a)  = r A.f  D^N  (a) 
A J J.k, 


Let  T f(x)  be  the  first  m terras  of  the  Taylor  series  for  f about  a, 
a 


i.e.. 


Tf(x)  = f(a)  + xDf(a)  + ...  + x 

a 


m-1  o'”  ^f(a) 


m-1  ! 


We  claim  that  X,f  = XT  f for  1 £ j £ m.  If  so,  then  (by  (3.13)) 

J J ^ 

D^F  f(a)  = D^F  T f(a),  0 < i < m-1.  Moreover  (Lemma  3.2),  since  T f is 
A A a — — a 

a polynomial  of  degree  m-1  < H,  F T f = T f.  Therefore 

A a a 

D^F  f(a)  = D^F  T f(a) 

A A a 


D^T  f(a) 
a 


D^f(a),  0 < i < m-l. 
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It  remains  only  to  verify  that  ^ T f for  1 < 1 < m. 

j j a _ J _ • 

According  to  (3.7)  and  (3.1) (a). 


'*'j(x) 


(a-x)  ^ ^Pj_i 


with  Pj  ^ a polynomial  of  degree  j-1.  Since  we  assume  that 


a. 


^(k-l-r) . 

J y 


4,(k-l-r)(a) 


0 for  all  j £ r < k. 


Therefore,  by  (3.6), 


=0  for  all  j < r < k, 
j,r  - 


and  ^jf  depends  only  on  the  first  j derivatives  of  f at  a.  But,  for 
1 < 1 < m,  the  value  of  these  derivatives  is  the  same  whether  T f or  f 
is  used. 


Corollary;  If  f(a)  = f(b)  =0,  = a,  and  = b,  then  F^f  ^ Sq(^,^). 

Later,  we  will  need  the  bound  on  the  weights  “ of  (3.6)  given  by 

J 

the  following  lemma. 

Lemma  3.4:  For  each  interval  I^,  1 £ i £ N,  for  all  j ^ V^, 

1“  1 £ C(k,M(^))  h^,  0 £ r < k, 

where  C(k,M(4))  depends  only  on  k and  M( 4) . 

Proof:  For  j e the  interval  (tj,tj^^^)  contains  I^.  Thus,  its 


PART  II 

CHAPTER  4 

PROBLEMS  WITH  0 < o < 1 


4.1  Introduction 

In  Chapters  4 - 8,  we  consider  the  singular  two-point  boundary 
value  problem 

(4.1)  - D(p(x)Du)  + q(x)u  = f(x),  0 < x < 1, 

(4.2)  u(0)  = u(l)  = 0, 

where  the  coefficients  p and  q satisfy 


a) 

p(x)  = x'’p(x) 

and  q(x)  = x°y(x),  where  0 £ o < 1, 

b) 

p(x)  > P . > 

0 for  all  X e I, 

~ min 

c) 

P e W^’"(I), 

d) 

Y e W°’°°(I). 
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Par  ter  [PI]  and  Greenspan  [G2]  considered  finite-difference  methods 
for  the  problem 

(4.4)  - D^u  + -iju  + qu  = 0, 

X 

u(0)  = 1,  u(l)  = 0. 

which  arises  by  separation  of  variables  in  the  equation  for  generalized 
axially  symmetric  potentials  in  a rectangle.  Problems  of  this  type  can 
be  replaced  by  equivalent  problems  of  the  form  (4.1)  - (4.2).  Jamet 
[Jl]  also  developed  finite-difference  methods  for  (4.4)  and  obtained 
L°°-norm  error  estimates  of  0(h^~°)  using  a uniform  mesh  of  size  h,  the 
exponent  l-o  being  sharp.  Gusman  and  Oganesyan  [G3]  considered 
finite-difference  methods  derived  from  variational  principles  for  the 
more  general  elliptic  problem 


( 

( 

Lu  = -1- 
3x 

3 

3y 

y^qCx.y) 

V 

I 

V 

in  a rectangular  domain  S = { (x,y)  | 0_<x£a,  0^y£b}  and  also 
obtained  low-order  error  estimates  using  a uniform  mesh.  Dershem  [D5] 
devised  second-order  accurate  three-point  finite-difference 
approximations  for  another  singular  ordinary  differential  operator, 

2 

Lu  = D ( xDu ) - — u . 

X 

Using  the  Rayleigh-Ritz-Galerkin  method  with  a singular  (L-spline) 
subspace  on  a uniform  mesh,  Ciarlet,  Natterer,  and  Varga  [Cl]  obtained 
L°“-error  estimates  of  order  h^~'’  (see  also  Chapter  8).  Natterer  [Nl] 
later  developed  a very  general  theory  for  approximations  of  this  type  by 
converting  the  problem  to  an  appropriate  singular  first-order  system. 


I 
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He  treated  problems  with  singularities  at  either  endpoint  and,  using  an 
appropriate  analogue  of  a B-graded  mesh  (with  grading  at  both  ends), 
obtained  Improved  error  bounds  of  order  N InCN)*-  with  c a positive 
constant  independent  of  N.  Other  subspaces  of  singular  functions  have 
been  considered  by  Crouzeix  and  Thomas  IC4]  (see  also  Chapter  8),  Jerome 
and  Pierce  [J2],  and  Dailey  and  Pierce  [Dl],  The  use  of  piecewise 
polynomial  spaces  on  a graded  mesh  was  also  suggested  by  Fried  and  Yang 
[F2]. 


4.2  Variational  Form  of  the  Problem 

In  this  section  we  instantiate  the  definitions  and  results  of 

Section  2.2  in  the  context  of  the  problem  (4.1)  - (4.2).  As  the 

0 2 

underlying  Hilbert  space  H,  we  use  H = L (1).  The  space  S,  in  which 
the  solution  of  (4.1)  - (4.2)  lies,  is  given  in  the  following 
definition. 

Definition  4.1:  Let  S be  the  linear  space  of  absolutely  continuous 
real-valued  functions  u on  1 such  that  u(0)  = u(l)  = 0 and 


Du  € (I ) . 


□ 


For  u,  V £ S,  we  define  the  inner  product 


a(u,v) 

and  the  norm 


1 

/ [p(x)DuDv  + q(x)uv]  dx 
0 


a(u,u) 


1/2 


J 
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1 

The  space  S is  an  example  of  a weighted  Sobolev  space  and  may  be 
identified  with  the  closure  of  the  C“  functions  with  compact  support  in 
1 with  respect  to  the  norm  ||  . 1|  if  the  problem  (4.1)  - (4.2)  has  a 
classical  solution  u,  then  u is  in  S [Cl],  but  is  not  necessarily  in 

H*^(I).  For  example,  let  p(x)  = x°,  q(x)  = 0,  and  f(x)  = 2-o.  Then  i 

u(x)  = X - X e H if  and  only  if  ° ^ "J* 


We  will  find  it  convenient  to  work  with  a different  norm  on  S, 


(4.5) 


S x°(Du)^  dx 

° J 


1/2 


u e S. 


We  assume  that  the  S-norm  and  energy  norm  are  equivalent,  i.e.,  that 


(4.6)  ^llvllg  < l|v||^  < A||v||g  for  all  V e S, 


where  ^ and  A positive  constants  Independent  of  v.  This  assumption  will 
be  satisfied  provided  q(x)  is  not  too  small,  in  particular,  if  -q(x)  is 
less  than  the  smallest  (positive)  eigenvalue  of  -D(pDu).  For  (a,b)  c 
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where  p “ (l-o) 

Proof : Let  v c S. 

|v(x) I 


and  * “ y* 


By  (A. 6), 


X 

< S |Dv(t) I dt 
0 


t'^/2 


dt 


< (/V'’  dt)l/2  (A^(Dv)2  dt)l/2 

0 0 

= Ullvllg 

£ <s  II V II 

E 

□ 


In  order  to  obtain  higher-order  error  bounds  for  the  solution  of 
(A.l)  - (4.2),  we  will  need  to  assume  that  the  right-hand  side  is  smooth 
(e.g.,  f e H™,  ra  ^0).  We  will  later  show  (Theorem  5.2)  that  D(x'’du)  is 
as  smooth  as  f and  therefore,  that  u Is  in  one  of  the  following  spaces. 


Definition  4. 2;  For  m ^ 0,  let 

s"  = {u  e S I D(x°Du)  e H®}, 

S ^ = {u  e S 1 x'^Du  e L“}, 

||D(x°Du)  II  = II  x°Du  II  ^00  for  u € S 

□ 

1-2 

The  role  of  the  space  in  Section  2.2  will  be  played  by  S 
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Let  u € S . By  Rolle's  theorem,  there  exists  x^  e i such  that 
DuCXq)  » 0.  Thus,  by  the  Cauchy-Schwarz  inequality. 


X Du(x)  = / D(t  Du)(t)  dt 

x„ 


£ 1|D(xDu)||q. 


Therefore,  for  u c S , 


(4.7)  l|D(x  Du)  II  II D (x  IIq  ^ •••i  IlDCxOi^) 


4.3  ^ B- graded  Meshes 

One  particular  kind  of  mesh  dominates  our  discussion  of  singular 
two-point  boundary  value  problems  in  Part  II.  These  are  the  "6-graded" 
meshes. 


Definition  4. 3:  Let  6^1  and  N a positive  integer.  The  6-graded  mesh 
Ag  is  the  partition  of  I given  by 


(4.8) 


0 < i < N. 


□ 


B-graded  meshes  were  introduced  by  Rice  [R2]  , who  considered  spline 
interpolation  of  the  function  x*,  o not  an  integer.  He  proved  that,  for 
6 appropriately  chosen  as  a function  of  a,  the  L“-error  in  interpolation 
on  a 6-graded  mesh  was  no  larger  than  that  obtained  when  Interpolating  a 
smooth  function,  despite  the  singularity  of  x*  at  0.  Eisenstat  and 
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mr 


Schultz  [E2]  have  applied  0-graded  meshes  to  a two-dljaenslonal  partial 

differential  equation  in  which  the  solution  has  a singularity  due  to  a 

corner  in  the  domain.  They  showed  that  convergence  at  the  best  possible 
2 

rate  in  L occurs,  despite  the  singularity,  when  using  a tensor  product 
mesh  which  is  graded  in  both  independent  variables.  Natterer  [Nl] 
employed  a 0-graded  mesh  for  a singular  two-point  boundary  value  problem 
for  a first  order  system  of  equations,  while  Fried  and  Yang  [F2] 
advocated  the  use  of  what  amounts  to  a B-graded  mesh  for  problems  like 
those  we  consider,  but  proved  no  error  bounds. 

For  B = 1,  the  mesh  A is  a uniform  mesh  with  h “ N As  B 

B ,N 

increases,  the  mesh  very  rapidly  becomes  skewed  towards  0;  to 
Illustrate,  we  depict  ^ in  Figure  4.1. 


i 

1 


LJ 1 ^ 1 I I 

0 X2  X3  X4  X5  1 

Figure  4.1:  A 

2,  6 

The  following  lemma  summarizes  several  important  properties  of 
B-graded  meshes. 
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k. 


Lemma  4.1:  There  exists  a constant  C(B)  = C(B,k),  independent  of  N, 

such  that  for  every  B > 1 and  all  N > 0,  A satisfies 

p ,N 

a)  h^  < Bi®"^N"®; 

b)  h < BN"^; 


(4.9) 


c)  M(Ae  < C(e); 


d)  le.l  < C(B)h.; 


e)  the  weights  id.  of  (3.6)  satisfy 

3 .r 


j,r' 


< C(B)h5^, 


for  all  j e V^. 


Proof:  By  the  mean  value  theorem,  there  exists  c e (i-l,i)  such  that 


. (1-1, » . 6c 


Thus , 


B B 

, 1 - (i-1) 

= 6 

N 


B-1  -B  „ B-1  -B 
< Be  N < Bi  " N , 


which  proves  (a).  By  the  same  argument, 
(4.10)  h^  > B(i-l)®"^  N"®  > h^_^. 


Thus,  the  sequence  (h^}  increases  monotonlcally  and 


h = h < BN®"^  N“®  = bn 

N - 


which 


-B 


proves  (b).  Clearly,  since  h^^  = x^  = N , 
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^2  2^-1^  6 
iT  = ~B—  “ 2 -1. 

"l  1 


while  for  2 < i £ N,  (4.10)  implies  that 


h,  B-1 

i < - " 


‘^1-1  (1-2)®  ^ 


< 3®-^ 


Therefore 


M(A  ) = max  = max(2®-l,  3®  ^) 

®’^  2<1<N  ‘^1-1 

(which  of  these  Is  larger  depends  on  B) , which  proves  (c) 


Next,  we  consider  I 9^^  I . For  all  1 £ i £ k, 


so  that 


'®ll  ^ ^i+k-1 


< (i+k-l)®h^ 


£ (2k-l)"h^. 


For  i > k,  9^  is  bounded  away  from  0,  i.e., 

®i  ’^i+k-l^  [Xpl], 


so  by  the  mean  value  theorem, 


'®1*  - ’'l+k-1  - ^1-k 


6 B 

= (i-Hk-l)  - (i-k) 

N® 


< (2k-l )B(i+k-l )®-^N-®. 
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CHAPTER  5 

PROPERTIES  OF  THE  SOLUTION 


5.1  Introduction 

In  this  chapter,  we  derive  several  properties  of  the  generalized 
solution  u of  (A.l)  - (4.2).  Our  goals  are  to  show  that 

1)  there  exists  a positive  constant  such  that,  if  u e S™,  then 

x“u  e and 

II  D^Cx^u)  II  Q < II  D*'"^  (x°Du)  II  Q,  1 _<  i < nrt-2; 

2)  there  exists  a positive  constant  F such  that,  if  f e H™  and  the 
coefficients  p and  y are  sufficiently  smooth,  then  u e S™  and 


||D(x'’Du)  11^  < r Ilf  11^; 

3)  if  u e S™,  then  there  exist  positive  constants  such  that, 
II  II  £ icj  II  D(x°Du)  II  -1  < i < m. 
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5.2  A Generalization  of  Hardy' s Inequality 

Let 

= {g  c C”(I)  I g(0)  = 0} 

and 

h”  = {g  € H^d)  I g(0)  = 0}. 

We  need  to  establish  several  properties  of  the  function  g/x  when  g e c” 

B 

or  g € ra  > 1.  In  Lemmas  5.1  and  5.3  we  show  that,  near  the  origin, 

D — 

g/x  behaves  essentially  like  Dg. 


Lemma  5.1:  If  g e m ^ 0,  then  g/x  e C™  and 

B 


(5.1) 

for  all  0 < 1 < m. 


x-»0+ 


Proof  I Clearly  g/x  e C^CO,!];  we  need  only  show  that,  for  0 ^ 1 £ m, 

11m  D^’Cg/x)  exists  and  has  the  value  given  by  (5.1).  For  x > 0, 
x-»0+ 

j, 

D*'(g/x)  - I (^)  D*’~j(x“l)  oJg 

j-0  J 

= E (J)  (-D^-Jd-j)!  oJg 

j=o  J 

= E If  (-l)^-J  xJ  oJg. 

j=0 


As  X 0,  the  last  sum  approaches  0.  Thus,  using  L'Hopltal's  rule. 
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lim  D*'(g/x) 
x-H)+ 


» lim 
x-»0+ 


I I 

Z ^ xJ  + E ^ jxJ“^  oJg 

j-0  J°1  

(1  + Dx*- 


= lim 
X+0+ 


iDl+lg  + E (-1)^-J+1  xJ-1  Djg+  E /.j.  (-D'-J  xJ- 

j=l  j=l 


=■  lim 
X-0+ 


D g 
4+r“ 


D^'^^g(O) 

~rfi — • 


□ 


In  order  to  extend  Lemma  5.1  to  the  case  g 


we  need  the 


following  result. 


Lemma  5.2:  For  all  m 1,  is  dense  in  Hg. 


Proof : Let  h e H™  c h™*  Given  e > 0,  we  can  (since  C™  is  dense  in  H™) 
B 

find  a g € such  that  II  h - g ||  < ■.  Furthermore,  by  the 

m l+L 

Sobolev  lemma  (Theorem  2.3), 

Cft 

llh  - gll^co  < ||h  - gll^,  < -f+cT"’ 


C,e 


1 . — 

since  m ^ 1.  Therefore  |g(0)|  < — . Let  g = g - g(0).  Then  g g C 

and 


8 - hll„  < Hi  - gll„  + llg  - h||„  < 


1 DJg 


J 


03  B 
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1 

Hardy's  Inequality  [HI]  states  that,  if  g e H^  then 

II  g/x  II  0 1 2 II  Dg  II  Q. 

The  following  lemma  generalizes  this  result  to  g e i 0. 

Lemma  5.3:  If  g e m ^ 0,  then  g/x  e H™  and 

(5.2)  IlD^g/x)  IIq  < ^ IlD^'^lgllg 

for  all  0 j<  i £ m. 

Proof:  We  demonstrate  (5.2)  for  g e The  result  then  follows  from 

" B 

the  denseness  of  in  H^^. 

We  first  show,  by  induction  on  I,  that 

(5.3)  xD^’Cg/x)  = D‘^g  - JlD''~^(g/x)  for  1 < 1 < m. 

The  case  J,  = 1 is  trivial.  Assuming  !,  > 1 , 

(5.4)  D(xD®'“^(g/x))  = xoVg/x)  + D*'“^(g/x), 

so  that 

xD*'(g/x)  = D (xD*'”^  (g/x) ) - D*'  ^(g/x). 

Using  the  induction  hypothesis, 

xD*'(g/x)  = DjD''"^g  - (4-l)D®-^(g/x)j  - D®'“^(g/x) 

= o'-g  - ID  (g/x). 

This  completes  the  induction. 

. J 
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Next,  note  that  for  1=0,  the  result  is  identical  to  Hardy's 
inequality.  For  1,  we  prove  (5.2)  by  computing  ||D*'(g/x)  1|  q. 
Integrating  by  parts, 

||D^(g/x)  II  2 , / x-2[xD‘'(g/x)]2  dx 

= -X  ^ [xD*'(g/x)]2  I 
'O 

+ 2 f x"MxD*'(g/x)]D  [xD®'(g/x)]  dx. 

0 

The  integrated  term  is  non-positive,  for 

- (D*'(g/x)  |^p2  < 0^ 

while 

lim  x(D*'(g/x)  = 0, 

x-*-0 

since,  by  Lemma  5.1,  D (g/x)  is  bounded  on  [0,1].  Thus 

l|D'^(g/x)  II  2 < 2 D''(g/x)  D(xD*'(g/x))  dx. 

0 

Using  the  identity  (5.3),  we  have 

lloNg/x)  ||2  < 2 D^g/x)  (D^+lg  - tD®'(g/x))  dx 

0 


= 2 / (D‘'(g/x)D^+lg)  dx  - 2M|D\g/x)||2 

0 ^ 

Hence,  by  the  Cauchy-Schwarz  inequality. 
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5 . 3 Regularity  of  the  Generalized  Solution 

We  turn  now  to  one  of  the  principal  results  of  this  chapter,  which 
will  later  be  used  to  show  that  solutions  of  (4.1)  - (4.2)  are  of  the 
form  x~*’v(x),  with  v a smooth  function.  This  result  is  crucial  to  the 
theory  of  Chapter  7,  in  which  we  consider  approximations  to  u of  the 
form  X with  v an  approximation  to  v.  Moreover,  it  provides  the 
means  for  proving  a regularity  result.  Theorem  5.2. 


Theorem  5.1:  If  u e S™,  m ^ 0,  then  x°u  e and  there  exists  a 

positive  constant  independent  of  u such  that 

(5.5)  ||d\x‘’u)||q  < llD^'^Cx^Du)  II  Q,  1 < i < mf2. 


Proof:  Define 


v.{t)  = t T>u(t)  , 0 < t < 1, 


(5.6) 


■^1+1(1)  = (v^(t)  - v^(0))/t,  0<t£l,  0^i<m. 


By  hypothesis,  Vq  e ; furthermore,  it  follows  from  Lemma  5.3  by  an 
obvious  inductive  argument  that  v.  e jjnri-l-i 


(5.7) 

provided  i < i £ m + 2. 


1-2 


ID^-Iy 


0 " 0 


We  now  show  by  induction  on  I that 


A,  o 


.i-l. 


i-2. 


(5.8)  D (x  u)  = D Vq  + c^D  + ...  + Cj^_^D 


''t-l(t)  dt. 


II 

11 

I: 

I 


L. 
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where  “ n (o  - i).  Since  u e S (and  in  particular,  u is 


absolutely  continuous) , 


x'^u(x)  = x°  / Du(t)  dt 


(5.9) 


=‘  x°  f t“°(t°Du(t))  dt. 


Differentiating  (5.9), 


T,/  0 x 0-1  f -0 

D(x  u)  = V + ox  f t V (t)  dt, 

u Q u 


which  is  (5.8)  for  1=1.  For  the  induction  step,  first  note  that  the 


last  term  on  the  right-hand  side  of  (5.8)  satisfies 


o-i  , 1-0 

X 0 


a-l  1-1-0  , 

= X ft  + tvj^(t)l  dt 


1 a-i  . ^l-a 

lU?''i-l(0)  + X ft  v^(t)  dt. 


Now  assume  (5.8)  holds  for  I,  Differentiating,  we  obtain 


„i+l,  a I i-1 

D (x  u)  = D Vq  + c^D  + ...  + Cj^_iDvj_^ 


, o-i-i  i-o 

X ft  V (t)  dt, 

0 


which  is  (5.8)  for  i+l. 


Inequality  (5.5)  now  follows  from  (5.7)  and  (5.8)  provided  we  can 


show  that 


J 


i-] 


"oil  7I-T^r 


which  is  just  what  we  need. 


o-l-  , , 
X I (x) 


0 II  ''t-1  II  0 


As  for  the  integrated  term,  at  x = 1 it  is  negative,  except 
(possibly)  when  !•  = 1 . But  when  1 = 1, 


I(x)  = / t'^VpiCt)  dt 

0 ^ 


/ Du(t)  dt, 
0 


whence 


1 

1(1)  = f Du(t)dt  * u(l)-u(0)  = 0; 

0 


(x)  dx. 
Then, 


so  the  integrated  term  vanishes. 
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At  X “ 0 it  vanishes,  too.  For  by  the  Cauchy— Schwarz  inequality. 


/ 


(I(x))^ 


/ t 


l-l-o 


''i-l(t)dt 


/ t2(^-l-°)  dt 

0 


J 

f 


\ 


^ vf  , dt 


V' 


i-1 


2i-20-l  X 

^ ir:2^ 


whence,  since  e L^(I), 

Lim  x2°"2M  (j- 
x-^0+ 


- 2l-2a-l  / '')i-l(t)  dt 

X-0+  0 


Q 


0. 


We  now  show  that  Theorem  5.1  applies  to  the  generalized  solution  of 
(4.1)  - (4.2)  provided  the  functions  f,  p,  and  y are  sufficiently 
smooth. 

Theorem  5.2;  Let  u e S be  the  generalized  solution  of  (4.1)-(4.2)  with 
f e H®,  m ^ 0.  Let  the  coefficients  p and  q satisfy  (4.3)  and,  for 
m > 0,  the  additional  hypotheses 


(5.11) 


P e ye  w®**. 


Then  u e S and  there  exists  a positive  constant  r.  Independent  of  f, 
such  that 


(5.12) 


D(x“du)  II  < r II  f IL,  0 < i < m. 


Proof:  We  start  by  showing  that  u c S®  and 
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II  D(x‘^Du)  II  < C II  f - qu  II  0 < t < m, 
where  C > 0 is  Independent  of  f. 

Following  Reddien  [Rl],  we  explicitly  construct  the  generalized 


solution. 

First,  let 

X 

(5.13) 

g(x)  = 

/ (qu- 
0 

■f)(t)  dt. 

Clearly  g 

e L“  and 

(5.14) 

II  g II  L- 

1 lU 

- qu  II  0* 

Next,  let 

(5.15) 

h(x)  = 

g(x)  = 
p(x) 

g(x) 

ptirr 

and 

X 

(5.16) 

w(x)  = 

f h(t) 
0 

dt. 

Finally, 

let 

(5.17) 

w(x)  = 

w(x)  - 

a(l)Y(x), 

where  Y 

is  the  solution 

of  the 

problem 

(5.18)  - D(pDY)  = 0 

(5.19)  Y(0)  = 0,  Y(l)  = 1, 

I 
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1 


We  claim  that  w e S and,  moreover,  w « u,  the  generalized  solution. 
First,  w(0)  = wQ)  = 0,  as  is  clear  from  (5.16),  (5.17),  and  (5.19). 
Second,  by  (5.16)  and  (5.17), 

Dw  = Dw  - w(l)DY  = 1 - w(l)DY, 

P 

so  that 

- C(1)x°/2dy. 

P 

Thus 

< -^||g|lj^„  llx"°^^  Hq  + lw(l)|  |1x°^^DY||q, 

min 

which  is  finite;  hence  w £ S. 


We  now  show  that  ||  u ” w 1|  ^ = 0,  and  hence  that  u = w.  Let  v e S. 
Integrating  by  parts,  using  (5.13)  and  (5.15)  - (5.18), 

1 1 
/ pDw(x)Dv(x)  dx  = / -D(pDw)v  dx 

0 0 

1 

= / [f-qu] (x)v(x)  dx 

0 

1 

= / pDu(x)Dv(x)  dx, 

0 

whence  taking  v = u - w , 

2 1 T 

l|u-w||p  = / p(D(u-w))^dx  = 0. 

^ 0 

This  shows  that  u = w. 
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By  (5.17), 


II  D(x  Du)  II  = II  D(x  Dw)  II 


< I|D(x'’dw)  II  + (0(1)1  IlDCxV) 


By  (5.15)  and  (5.16), 


d'-'^xV)  = D‘+^x"h)  = d'+NI).  0<l<m. 


.^+1/.A 


whence 


II  D(x  Dw)  II  < 


I 1 CQ  II  g II 


£ Cj^  II  f - qu  II  0 £ *•  £ m, 

where  Cj^  depends  only  on  and  l|P||y®',°°.  Also,  by  (5.1A)  - (5.16), 


I w(l ) I £ / I h( t) I dt 

0 


g II  , 00  / x"  dt 


- (l-o)P  " ® " L” 
min 


- (T-a-frr  " ^ " 0- 


||D(x°Du)||^  < ||D(x‘’dO)  II  + |w(l)|  ||D(xV)  II 


£ C II  f - qu  II 


where  C » inax(cp  - 
independent  of  u.) 


||D(x  DY)  II  ^ 


).  (Note  that  ||D(x  DY)  ||  , is  a constant. 


r 
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Now  we  show  that  there  exist  constants  such  that 
(5.21)  Ilf  - quit  ^ Ilf  II 

We  proceed  by  induction  on  For  i » 0,  the  a priori  estimate  of  (2.7) 
yields 

lUllg  < Ilf  No, 

whence 

II  f - qu  II  Q < (1  + 6^  llq  II  II  f II  ^ 

= «o  llfllo- 

Since  DCx'^Du)  e we  can  apply  Theorem  5.1  which  yields 
u = X °v,  V e 

and 

II V 11 2 < ||D(xV)||q  < F^cMq  IlfllQ. 

Thus  qu  = (x°y) (x”°v)  “ yv  and 

II  f - qu  II  2 < II  f II  2 + II  Yv  II  2 
1 Ilf  II  2 + ^2  II  V II  2 

< ”2  Ilf  II  2’ 

where  C2  depends  on  II  Y 11^2,”  and  M2  - 1 + C2fj^C^MQ). 

We  reiterate  this  argument  for  1-2,4,  ...,  m (if  m is  even)  or 
1-2,4,  ...,  m-1  (if  m is  odd).  To  obtain  (5.21)  for  odd  values  of  t, 
possibly  including  m,  we  apply  Theorem  5.1.  Since  D(x^Du)  € 

X V,  where  v £ H , and 


J 


u 


^1 

i 


□ 

5. 4 Bounds  on  Weighted  Norms  of  the  Generalized  Solution 

The  last  two  results  of  this  chapter  establish  bounds  on  two 
"weighted  norms"  of  the  derivatives  of  the  solution. 

Lemma  5.4:  There  exist  constants  which  depend  only  on  I and  o such 
that.  If  u e s“,  m ^ -1 , then  x°'^*D*''''^u  e L“(I)  and 

(5.22)  11^.  < ||d(x‘’Du)  II  j 


for  all  -I  < !.  < m. 
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I V II  1+1  < IlDCx^Du)  II 


Thus 


1 r c II  f - qu 


1-1 


i ^^1-1  iinit-i. 


f - qu  II 1 < II  f 11 1 + 11  Yv 


< II  f II  j,  + Cl  II  V II 1 


1 II  f II  1 + Cl  II  V II  1+1 


1 M,  Ilf  11^, 


where  c^  depends  on  ||  y 11^8-,“  and  Mi  = 1 + CiriCMj_l). 
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Proof;  First  take  t - -1.  By  the  definition  of  S, 
iu(x)i  < iDu(t) I dt  - -tUHnilLL 


dt 


0 t 


< II  x°Du  II 


Therefore 


|x°”^u(x)|  < -pL-  II  x°Du 


IlDCxV)  ||_^. 


For  1^0,  we  use  induction  on  i.  Since 

t+l  a . Or.*'+l  o-i 

D (x  u)  = X D u + I *^4  X ^ D Ju, 

j-1  ^ 

1+1 

where  c = ( ^^ ) n (a-i),  it  follows  from  the  induction  hypothesis 
•J  J 1=0 

that 


1+1 


(5.23) 


l^a+lpl+i^l  ^ |x*'D*'+^x‘’u)  I + Z |c  D^’+^'Jul 

j=l  ^ 


.1+1  . <J 


1+1 


< l|D  (xu)||^=+  Z Icjl  II  D(x  Du)  II 

By  Sobolev's  lemma  (Theorem  2.3), 


II  D^‘^\x‘'u)  II  - ||d\d(x'’u))  llj^a 


< CJ|D(x  u)  11^^^, 


and  by  Theorem  5.1, 


Thus 


II  D(x°u)  II  1^1  < Ti  II  D(x‘’Du)  II 


II  D*’’*'*^  (x^u)  II  _<  II  D(x‘^Du)  II 
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Together  with  (4.7)  and  (5.23),  this  yields  (5.22)  with 

t+1 

Z 

J = 1 


G 


Lemma  5.5:  For  each  positive  integer  1 and  every  e > 0,  there  exists  a 
constant  = K(o,e,J.).  independent  of  u,  such  that,  if  u € S®,  m ^-1, 


then 


(5.24)  l|x(2^-3+2‘’+^)/2  D^ll  < K Ij  D(x‘’Du)  1L_ 


1-2 


for  all  1 < 1 < nri-2. 


Proof:  Again  we  use  induction  on  i.  For  <■  = 1 , 


1 


1 


^“^^(00)^  dx  < ||x°Du  II  2.  dx 

0 ^0 


= e“Ml  D(x°Du)  " 2 


r 


For  I > 1,  we  differentiate  to  obtain 


t-1 


whence 


Let 


d''-1(x°Du)  = x^dV  + ”'r\'':SD^(x'’)D^-iu, 

i=l  ^ 


x'^D^'u  = d’^"2[d(x‘’du)]  - (^7^)D^(x°)D*'"^u. 

i=l  ^ 


= ( .S  n (0-j), 

j=0 


r 
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By  the  triangle  inequality, 

||^(2»-3t2<.+c)/2  11^ 

, |U<2‘-3«)/2 

t IcJ  ||,<2(t-l)-3+2»f.)/2 

i-1  ^ 

Since  i ^ 2,  2J.  - 3 + e > 0,  so  that 

||x(2‘'‘^+^)/^  d'^"2(D(x'’Du))  II  Q < ||D*’"^(D(x‘’Du))  II  Q. 

Thus,  by  the  inductive  hypothesis, 

„^(2i-3+2o+e)/2  ,|  ^ 


i-1 

£ (1  + Z |c  |K  ) IlDCx^Du)  II  2» 


CHAPTER  6 


SPLINES  ON  A 6-GRADED  MESH 


6. 1 Introduction 


In  this  chapter  we  explore  the  use  of  splines  on  B-graded  meshes  to 
approximate  the  solution  of  (4.1)  - (4.2).  Subspaces  of  piecewise 
polynomials  have  been  studied  extensively,  and  their  properties  as 
approximators  of  smooth  functions  are  well-known  (see,  e.g.,  [S3], 

[D2]).  Our  aim  here  Is  to  find  analogues  of  these  results  for  the  case 
of  functions,  elements  of  s'”,  with  a singularity  at  0. 

Given  u e s'”,  we  wish  to  show  that  with  an  appropriately  chosen 
mesh  the  quasllnterpolant  Is  a "good"  approximation  to  u,  even  though 
u ^ H™  for  any  m ^ 2.  We  cope  with  the  singularity  In  u by  grading  the 
mesh  so  that  the  Intervals  near  zero  are  small.  We  employ  a B-graded 
mesh,  the  definition  and  properties  of  which  are  given  In  Section  4.3. 


In  Section  6.2  we  prove  error  bounds  for  the  quasllnterpolant  of 
functions  In  s'”.  S-norm  bounds  are  obtained  In  Theorem  6.1  and  L^-norm 
bounds  In  Theorem  6.2.  In  both  theorems,  we  specify  a B-gradlng  which 
yields  optimal-order  accuracy.  Let  k ^ 2 be  the  order  of  the 
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We  obtain  the  following  optimal-order  error  bounds  for  the 
quasllnterpolant  F^u  e 

llu-F^ullg  < C ||D(x°Du)  II  ^_2  ife>Bi. 

II  u - F^u  II  Q < C N II  DCx^Du)  II  if  B > 03. 

Obviously,  since  I > 2,  6^  < 82  < 6^.  Moreover,  for  fixed  i,  6^  and 

B2  ® as  o 1,  while  B3  remains  bounded  for  all  0 ^ 0 < 1.  Thus,  far 

2 

less  B”grading  is  required  to  prove  the  L -error  bound. 

As  a corollary  to  these  approximation-theoretic  results  we  obtain 
2 

(in  Theorem  6.3)  L -norm  and  S-norm  error  bounds  for  the  RRG 

approximation  to  the  generalized  solution  of  (4.1)  - (4.2),  provided  the 

heavy  B'grading  of  Theorem  6.1  (b  > Bj^)  is  used.  In  Section  6.3  we  show 

that  the  RRG  approximation  will  not  achieve  the  optimal  rate  of 

2 2 

convergence  in  the  L -norm  in  the  spaces  S (a  ).  We  conjecture  that 

U Bo,N 

■'  2 

while  B “ Bg  is  too  small,  B “ 62  suffices  for  optimal  L -norm  accuracy 
of  the  RRG  approximation,  while  B “ B^  is  required  for  optimal  S-norra 
accuracy.  The  numerical  results  in  Section  6.4  Include  an  empirical 
study  of  the  effect  of  the  grading  parameter  b on  the  error  and  rate  of 


convergence  of  the  RRG  approximation  which  tends  to  substantiate  this 
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conjecture.  Other  experiments  confirm  the  approximation-theoretic 
results  of  Section  6.2.  In  Section  6.5  we  consider  questions  of 
computational  complexity. 

6.2  The  Quasiinterpolant  in 

j Our  proof  that  splines  can  approximate  a function  u(x)  which 

i behaves  like  x^  ° depends  on  the  existence  of  local  approximation 

I mappings  (for  example,  the  quasiinterpolant).  For  such  an  approximation 

mapping,  the  error  in  the  interval  depends  only  on  the  behavior  of  u 
I in  a small  number  of  neighboring  intervals.  Thus,  by  using  a B-graded 

i mesh,  we  overcome  the  growth  of  the  derivatives  of  u near  0. 

I 

We  want  the  quasiinterpolant  to  satisfy  the  boundary  conditions  of 
(4.2).  By  the  Corollary  to  Lemma  3.3,  if  we  choose 

I (6.1a)  ''l  “ ° ^d  ° 

then  F^u  e Sq(A,£)  for  all  u e S.  Here,  and  throughout  the  rest  of  the 
chapter,  ^ may  be  any  Incidence  vector. 

We  would  like  to  make  the  quasiinterpolant  independent  of  u near  0 
for  all  but  the  first  subinterval.  Therefore,  we  require  that 

' (6.1b)  i > 2 ^ j _<  d. 

I . 

I 

Thus,  in  the  intervals  I^,  2 < i £ N,  6^  c [x^,l].  The  error  bounds  of 


Theorem  3.1  then  apply,  since  u is  smooth  in  0^. 
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In  the  first  interval.  Theorem  3.1  is  of  no  help,  since 

I 2 

D u ^ L (0,x^)  for  any  i 1.  But  since  the  mesh  is  6-graded,  so  that 
a 

Xj^  - N~  is  very  small,  it  is  possible  to  bound  ||  u - ^^u||g(jj  by 

bounding  first  11^*11  s(l)»  then  llf^ullg(j),  and  using  the  triangle 
inequality. 

We  first  obtain  bounds  on  the  linear  functionals  ^jU,  for  j c V^. 


Lemma  6.1;  There  exists  a positive  constant  Cj^  Independent  of  u and  N 
1-2 

such  that,  if  u e S , I ^ 2,  A = and  the  quasiinterpolation 

points  satisfy  (6.1),  then 

(6.2)  1 ‘'I  IID(x''Du)  II 

for  all  2 ^ j £ k. 

Proof ; Using  (6.1b)  and  the  bound  of  (4.9)(e), 

l^.ul  = z [0.4  r d’^u(_)  I 

-*  r<!. 

< C(B)  z hi  |d>^u(  )| 

r<i 

< C(S)  x}-'  j 1 Dr„(  j_ 

r<i 

since  h^  “ *1*  Now,  using  Lemma  5.4  and  (4.7), 


^.u|  _<  C(B)  x,“  Z xr_i  II  D(x  Du)  II 

J r<l 

1 C(6)  ( z «__f)  x}"'’  IID(x°Du)  \\^_2' 
r<t 
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Using  this  result,  we  obtain  S-nonn  and  L -norm  bounds  on  F^u  in 
the  first  interval. 


Lemma  6.2;  There  exists  a positive  constant  Independent  of  u and  N 
1-2 

such  that,  if  u e S . ^ i 2,  A = and  the  quasiinterpolation 

points  ^■'j}  satisfy  (6.1),  then  F^u  = 2 t ^ satisfies 


(6.3) 


(6  4) 


1 i l|D(x°Du)  II  ,_2. 


Proof:  In  the  proof  of  Lemma  3.3,  it  is  shown  that  X^u  » u(0)  =•  0. 
Therefore,  in  Ij^,  F^u  is  given  by 


Thus , 


'a'^  II  S(l^)  1 l^j^l  ll’^j.k  II  S(Ip 

k. 

'A'^llL'd^)  - ^^2  llNj,k  IlLdl^)* 


Using  the  result  of  Lemma  3.1, 


■ j,k  " S(I^) 


f x'’(DNj^j^)^  dx 


< II DN 


j,k  "L"(ip  1+5“ 


1 1/2 


(1+0) 


— — “ Ai 
1/2  1 


while 


1 
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/ Xi 


V 

< II  N 


/ (N  dx 

0 ^ ’ 


1/2 


< . 


j.k  "l*(I^)  ’'I 

1/2 


1/2 


Thus,  by  (6.2), 


„(l-o)/2 

(1+a)^ 


1 7^7:717^ 


1-2 


and 


- (k-l)C]^  x^^^  II  D(x“Du)  II  j^_2. 


We  now  obtain  S-norm  error  bounds  for  the  quasiinterpolant. 


Theorem  6. 1;  If  ue  S™,m^0,  A = AgjjWith6>  6^  = 


2(t-l) 


evaluation  points  i’^j}  satisfy  (6.1),  then  F^u  “ 2 i ^ Sq(A,^)  and 

(6.5) 


II  S - C3N~^*"^^  II  D(x'’du)  II  j^_2, 


where  I = mln(k,  iiri-2 ) and  is  independent  of  u and  N. 


Proof : We  consider  Ij^  separately  from  the  other  intervals.  By 
Lemma  5.4, 

X,  ^ ^ 

x^^Dul^dx  = r iiL 

0 


u||^,T  N = dx 

0 0 


1-<J 


— ^ CT 


I x'^Du  II  ^00 


< x|"°  ||D(x°Du)  11^. 


J 
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By  (4.8),  xp  ®)/2  _ and  by  the  triangle 

inequality  and  (6.3), 


(6.6) 


Hu 


4^I‘S(I^)  1 


< 


I*  “ I'  s(ip  II  sdp 

(U<Q  + C2)  ||D(x°Du)  II  ^_2 

+ C2)  II  D(x°Du)  II  j,_2 

||D(x‘'Du)  II  j_2. 


Now  we  consider  the  remaining  intervals  I^,  2 < 1 < N.  The  error 
bounds  of  Theorem  3.1  refer  to  the  neighborhood  9^  of  1^.  Let  be  the 

largest  integer  such  that  0^  c [x.  ,1].  By  the  definition  of  9.  and  the 
conditions  (6.1)  on  the  points 


j ^ max(l,  1-k) 


(see  Figure  6.1). 


Figure  6.1:  Definition  of 
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and  j-  < k,  since  > maxd,  i-k).  Thus 


u - F 


4“  II  S(I^) 


(6.7) 


< K ; x'^(D^u)2  dx. 

e. 


c2  t^-2(i-l)  ^ x^(D^u)2  dx. 


The  intervals  9^  overlap,  but  no  more  than  2k  of  them  at  any  point. 
Thus,  by  Lemma  5.5, 


||x^/2Di^  ^ 2k 

< 2k  k2  ||d(x°Du)  ||2_2. 


(By  the  assumptions  on  6,  e > 0 as  required  by  Lemma  5.5.)  Together 
with  (6.6)  and  (6.7),  this  yields 


r ti  ^ 
II  S 


II"  -^au1Is(i^) 


< (c2  + 2k(K^C5)^)  N-2(‘-n  H D(x“Du)  |1  2_2. 


We  now  consider  L -norm  error  bounds. 


Theorem  If  u e S , m > 0,  A = A with  6 > e,  = ■ v , and  the 

B ,N  3 3 - 2a 

evaluation  points  {t^}  satisfy  (6.1),  then  F^u  = F^^^  t ^ So(A,z) 


and 
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(6.8) 


“ II  D(x°Du)  II  j^_2. 


where  I = min(k,  nH-2)  and  is  Independent  of  u and  N. 

Proof : Again  we  consider  separately  from  the  other  intervals.  By 

(4.8),  « f}-®(3“2°)  < n~2*'^  Lemma  5.4  and  (4.7), 

- 

II  u II  ^2  » ; dx 

L up  0 


/ (x'^'^^u)^  dx 


(6.9) 


< ^ x3-2°  ||x-lu||2. 


< 3^  II  D(x°Du)  II 


- 1=2S 


JJ  2 


By  the  triangle  inequality,  (6.4),  and  (6.9), 


(6.10)  ll'^  - V 1 (C2  + -^)  N~^  ||D(x°Du) 

^ />2? 

= Cy  n"*'  I|D(x®Du)  II  j^_2. 


We  now  consider  I^,  2 ^ 1 £ N.  According  to  Theorem  3.1  and 


(4.9)(d), 


/ (u  - jjx  < k2  le^l^*"  / (D^-ul^  dx 


< 0(6)^*'  h^*'  / (D*'u)^  dx. 

®i 
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Because  B(3  - 2a)  > 2l, 


e E (3  - 2o)  - 21/6  > 0 


6(3  - 2o  - e)  = 21. 


Let  X and  C be  defined  as  In  the  proof  of  Theorem  6.1.  Then 


; (u  - F u)^  dx  < 0(6)^“'  ~ f x'^(d‘'u)2  dx 


(6.11) 


< (6C(B))^*'  ; x^(D®'u)^dx 


< (6C(6)k®“M^*'  N"^*'  ; x'^(D*'u)^  dx 


c / x^(D®'u)2  dx, 


since  1/j^  < k,  B(2«-t)  = 6(3  _ 2a  - e)  = 2«-,  and 


6C  = 6(21  - 3 + 2o  + e)  = 216  - 6(3  - 2o  - e)  = 21(6  - 1). 


From  (6.10)  and  (6.11), 


N ’'i 


/ (u  - F dx  = I f (u  - F dx 


i=l  X, 


£ N c^llD(x'^Du)  11  0 , + Co  I J x^(D*’u)^  dx  . 

I ‘-2  ^ 7 


more  than  2k  of  the  0^  overlap  at  any  point,  so  by  Lemma  5.5, 


1 

! 

i 
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Z ; x^(D^u)2  dx  < 2k  ;^x'^(d‘u)2  dx 

1-2  0.  0 

2 I 

< 2kK^_2  ||D(x°Du)  ||2_2.  \ 

I 

This  yields  (6.8),  with  - (c^  + 2kKj_2C8 

Q 

We  now  obtain  error  bounds  for  the  Rayleigh-RltzMJalerkln  procedure 
for  approximating  the  solution  of  the  singular  problem  (4.1)  - (4.2). 

Theorem  6. 3:  Let  u c approximation  to  the 

generalized  solution  u of  (4.1)  - (4.2),  where  6 > B = — , if  p 

1 1-0 

q,  and  f satisfy  the  hypotheses  of  Theorem  5.2,  then 

(6.12)  Hu  - till  g < X ‘^3r  Ilf  llt-2 

and 

(6.13)  Hu  - till  0 < (Ac3r)2  N"*-  H f 11  j^_2, 

where  I = min(k,  nH-2). 

Proof:  By  Theorem  6.1  (the  "approximation  hypothesis").  Theorem  5.2 
(the  "regularity  hypothesis"),  and  Nltsche's  trick  (Theorem  2.6), 

II  “ ■ “ II  S - X *^3  H D(x®Du)  II  j_2 

and 

Hu  - till  Q < (Ac3)^  r N“^  HD(x‘^Du)  ||j_2. 

We  now  use  Theorem  5.2  to  bound  the  right-hand  side  of  these 
Inequalities. 


□ 
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6.2  A Counterexample 

To  the  best  of  our  knowledge,  no  example  has  yet  been  reported  In 

which  the  Rayleigh-Rltz-Galerkin  method  failed  to  produce  an 
2 

optimal-order  L -norm  approximation.  Eisenstat,  Schultz,  and  the  author 

[El]  have  shown  that  no  such  example  exists  for  nonsingular  problems. 

For  the  differential  equation  (4.1)  - (4.2)  in  which  p > p >0  and  a 

~ min 

family  of  C^-piecewise  polynomial  subspaces  with  respect  to  a 

quasiunlform  mesh. 


0 i 


inf 


^h  " 0- 


where  u is  the  RRG  approximation  to  u in  S,  and  the  constant  C is 

h 

independent  of  u and  h.  However,  for  the  differential  equation  which 

2 

concerns  us  here,  there  do  exist  subspaces  for  which  the  L -norm  of  the 

error  of  the  RRG  approximations  do  not  decrease  as  fast,  as  a function 

2 

of  h,  as  the  error  in  the  best  L -approximations.  We  shall  now  give 
such  an  example. 


We  restrict  our  attention  to  the  simple  case  of  continuous, 

2 

piecewise  linear  splines  (the  subspace  S (4),  which  we  hereafter  refer 

to  as  L(4)).  Let  u be  the  generalized  solution  of  (4.1)  - (4.2). 

According  to  Theorem  6.3,  if  D(x°Du)  e and  u e L(4  ) ig  the  RRG 

p ,N 

2 

approximation  to  u,  where  6 > 6^  = yj,  then 

II  u - u II  Q < CN  ^ II  D(x  Du)  II  Q. 

This  much  6-grading  is  not  required  for  a 2”‘^-order  L^-norm 
approximation  to  exist.  The  result  of  Theorem  6.2  is  the  error  bound 
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llu  - UiIIq  < CN”^  llDCx^Du)  II  Q, 
where  e L(^g  jj)  is  the  quasllnterpolant  of  u and  6 > - '^~2a* 

If  we  choose  as  the  quasiinterpolation  points  of  (3.5)  - (3.6) 

(6-lA)  Tj  - Xj_^,  1 < j < d - N+1, 


then  the  quasllnterpolant  Interpolates: 


^1^*1^  ” u(xj^) , 0 ^ 1 £ N. 

In  other  words,  the  piecewise  linear  quasllnterpolant  with  the  "obvious" 
choice  (6.14)  for  Tj  is  the  familiar  piecewise  linear  interpolate.  To 
show  this,  we  note  that  for  k » 2,  the  quasiinterpolation  formulae 
involve  only  u(Tj)  and  Du(ij).  By  choosing  Tj  = Xj_j^,  we  make  the 
coefficients  “j  “ 0 and  q = 1;  thus  “ u(Xj_j^). 

Let  us  define  the  "hat  functions"  c L(4)  by 


Li(xj)  = 0 < i,j  < N 

(see  Figure  6.2). 


^1-1  *i  ^i+1 

Figure  6.2:  The  "hat  function"  L^(x). 


It  is  easily  verified  that  2»  the  "hat  functions"  are 

the  piecewise  linear  B-spllnes.  We  have  shown,  therefore,  that  the 
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quaslinterpolant  is  given  by 

N+1 


and  that 


N N-1 

= I u(x  )L  = I u(x  )L 
1=0  1=1 


Uj(x^)  = u(x^)Lj(x^)  = u(x^). 


In  order  to  show  that  the  RRG  approximation  is  not  optimal  in  this 
case,  we  need  a lower  bound  on  the  error  Hu  - u ||  q.  But  since 

II  U - U II  Q ^ II  U “ Uj  H Q “ II  U - Uj  II  Q 

and  II  u - Uj  II  Q =0(n“^).  it  suffices  to  prove  that  ||  u - ||  q > 

for  some  C,  e > 0.  We  will  in  fact  show  that  ||  u - II  q > CN”^» 

Our  plan  of  attack  is  as  follows.  We  know  that  the  coefficient 
vector  of  u with  respect  to  the  basis  {L^}  is  the  solution  of  the  linear 
system  A^=  f^  of  (2.10).  For  a specific  problem  (that  is,  a specific 
function  f and  corresponding  solution  u) , we  will  compute  the  elements 
of  A and  (Fortunately,  A is  tridiagonal.)  Since  we  know  u,  we  also 
can  compute  the  coefficient  vector  of  u^  and  the  vector  £ = Ah^.  It 
happens  that  the  elements  of  ^ ^ are  positive  and  that  A ^ has 

positive  elements.  Thus  the  elements  ofw  = ^-  _n=A  ^(^-_g)  are  also 

positive.  We  then  get  a lower  bound  on  the  norm  of  u - u by  exhibiting 

2 

a vector  _C  with  the  property  that  A^  < Aw  and  showing  that  the  L -norm 
of  the  element  of  L(A)  with  coefficient  vector  t decreases  at  the  rate 


0(N"^) 
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r 
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dx  = / dx 

’‘i-1  ’'i-1 


(6. 16) 


= 1 j^-2  3/2  _ 3/2 

3 i ^ i i-r 


= ^Ni 


^ - (i-1)^ 


and 


3 2 

(21-1)^ 


3 ^ 3 

-2  ^ L^(x)dx  = (h^  + 


3iN 


From  (6.15)  and  (6.16), 

Ail  = £ = . 


where 


8. 


■‘’iVi  + + - ^+iVi 


3iN"^  1 - 


1/3 


3(i6i^  - 8i^  + 


4i 


3(16i^  - 8i^  + 1) 


N-1 

N-1 

Let  e “ u - u = ^ (^. 

^ i=l 

- ^,)L^  = ^ 

^ i=l 

Vi* 

Then 

Aw  = £ “ f.  = 

4i  - iN“^ 

— > 

n 

3(16i^  - 8i^  + 

1) 

yj  • 

In  order  to  establish  a 

C(N“h  + o(N”h 

lower 

bound 

N-1 


employ  the  piecewise  linear  function  e = f ^i^i  ^ ^0^^^  given  by 

i=l 
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Cl  = —h  1<1<N-1. 

We  claim  that  II  e ||  ^ - CN"^  +0(N"^/^)  and 

(6.17)  0 < e(x^)  < e(x^),  1 < i < N-1 

(i.e.,  0 < < w^,  1 £ i £ N-1). 

In  order  to  show  (6.17),  we  make  use  of  the  fact  that  A is  an 
irreducible  Stieltjes  matrix  Thus  [VI,  Corollary  3 of  Theorem  3.11), 
A”^  > 0.  This  implies  that,  for  any  n-vectors  ^ and  £,  if  A^  > An^,  then 
£ > _n. 

First  we  compute  II  e ||  q.  Note  that,  in  the  interval  [xpl],  e is 
monotonically  decreasing.  Thus 

lie  11^  > 2 h.  e(x.)^ 

^ i=2 

i=2  N^  ' 6N^  ' 

1 N 2 

= I (2i-l)(N-i)^ 

36n'’  i=2 


A Stieltjes  matrix  is  a symmetric  positive  definite  matrix  with 
non-positive  off-diagonal  elements.  As  in  [VI],  we  say  for  n x n 
matrices  A = [^^j]  » B = that  A > B if  a^^j  > b^j  for  all 

1 _<  i,j  £ n.  We  say  A > 0 if  a^^  > 0 for  all  1 _<  i,j  £ n.  These 

notational  definitions  also  apply  to  n-vectors,  considered  as  n x 1 
matrices. 
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As  we  noted  above,  in  order  to  prove  (6.17)  we  need  to  compute 
and  check  that  A^  < Aw  = £ - f.* 

(A^)j^  * (bj^  + h2  ) ^ “ ^2^2 


-2 

since  (Aw)^^  ..  ^ .1. 


Fo-  2 _<  i £ N-1,  < b^,  whence 

(AO^  - [-b.(N-i+l)  + (b^  + b.^^)(N-l)  - b^^^(N-i-l)] 

6N^ 

= - bi)  < 0 < (Aw)^. 

We  have  shown  that  II  ® II  p > II  ® II  q “ CN  ^ +0(N 

Therefore, 

I|u-u||q  > Ilu-Ujll^-  Ilu-Ujilo  > CN“^  + 0(N-3/2)^ 

Thus,  in  this  situation,  the  RRG  procedure  does  not  produce  an  optimal 
2 

L -approximation. 
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6.4  Numerical  Results 

In  this  section  we  present  the  results  of  some  numerical 
experiments  which  illustrate  the  computational  procedure  analyzed  in  the 
previous  sections.  We  consider  the  problem 


- D(/x  Du)  = 


TT^  sin(7rx)  + -^(  it  cos(irx)  - ), 


u(0)  = u(l)  = 0, 


which  has  the  solution  u(x)  = ^ approximations  to  u 

from  various  spline  subspaces  were  computed  and  the  error  tabulated 
below. 

All  computations  were  performed  in  double  precision  on  a PDP-10 
(with  54  binary  digits).  The  integrals  required  were  computed  using 
Gaussian  quadrature,  with  k-1  nodes  in  each  interval  of  the  mesh.  In 

the  tables  below,  e^^  denotes  the  error  u - u,  where  u is  the  RRG 

k 2 

approximation  in  j^).  We  measure  the  error  in  the  L -norm 

( II  II  q)  S-norm  ( ||ejjl|g),  which  were  computed  using  Gaussian 

quadrature  with  k+1  nodes  in  each  interval.  We  also  computed  the 

quantity 


'N  " “,4 


“ax  Ie  (x)l. 
l<i<N 


In  the  following  tables,  we  use  the  shorthand  notation 


-4 


5.2  (-4)  for  5.2  x 10 


Next  to  the  error  e for  a given  value  of  N,  we  list  the  observed  rate 
N 
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N lleNllo  rate  lleNll-.A  RATE  II  II  s ^‘ATE 


16 

.79 

(-02) 

.11 

(-01) 

.12 

(+00) 

32 

.89 

(-03) 

3. 15 

.11 

(-02) 

3.32 

.30 

(-01) 

1.98 

48 

.25 

(-03) 

3. 16 

.22 

(-03) 

4.02 

.13 

(-01) 

2.00 

64 

.10 

(-03) 

3.09 

.64 

(-04) 

4.22 

.75 

(-02) 

1.98 

80 

.52 

(-04) 

3.06 

.25 

(-04) 

4.29 

.48 

(-02) 

1.98 

96 

.30 

(-04) 

3.04 

.11 

(-04) 

4.31 

.34 

(-02) 

1.97 

Table  6.2: 

Error 

in  Sq(A 

8,N^  ® 

= ®1 

N 

II 

II  0 

RATE 

II 

II  m,  ^ 

RATE 

II  eN 

II  S 

RATE 

16 

.97 

— 

(-03) 

.11 

(-02) 

.23 

(-01) 

32 

.12 

(-03) 

3.01 

.20 

(-03) 

2.36 

.40 

(-02) 

2.50 

48 

.31 

(-04) 

3.35 

.55 

(-04) 

3.22 

.14 

(-02) 

2.59 

64 

.91 

(-05) 

4.27 

.15 

(-04) 

4.44 

.63 

(-03) 

2.81 

80 

.34 

(-05) 

4.45 

.54 

(-05) 

4.75 

.33 

(-03) 

2.85 

96 

.15 

(-05) 

4.37 

.34 

(-05) 

2.42 

.20 

(-03) 

2.85 

Table  6.3:  Error  in  jj)  — 12  » 
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In  order  to  get  some  idea  of  the  effect  of  the  parameter  6 on  the 
error  H II  q,  we  computed  the  RRG  approximations  in  several  cubic 
spline  subspaces  (the  spaces  Sq(A)).  we  used  meshes  with  N - 16,  32, 
and  64  and  6=1,  2,  ...,  13.  These  results  are  presented  in  Figure 

6.3. 


The  results  of  this  experiment  lead  us  to  believe  that  the  RRG 

Ic  ^ ^ 

approximation  in  S (A  ) converges  in  the  L^-norm  at  the  rate  N“*^  and 

that  the  error  will  generally  be  smaller  than  for  the  more  heavily 


graded  meshes  with  B = The  results  of  an  experiment  to  test  this 

hypothesis  are  presented  in  Tables  6.4  - 6.6. 
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Table  6.4:  Error  In  S 


.53  (-04) 
.79  (-05) 
.36  (-05) 
.19  (-05) 
.11  (-05) 


32 

.38 

— 

(-03) 

3. 18 

.36 

— 

(-03) 

— 

3.79 

.21 

— 

(-01) 

1.85 

48 

.11 

(-03) 

3.08 

.73 

(-04) 

3.95 

.10 

(-01) 

1.76 

64 

.46 

(-04) 

3.04 

.32 

(-04) 

2.87 

.61 

(-02) 

1.71 

80 

.23 

(-04) 

3.02 

.17 

(-04) 

2.90 

.42 

(-02) 

1.68 

96 

.13 

(-04) 

3.01 

.99 

(-05) 

2.91 

.31 

(-02) 

1.65 

— 6 = 6, 


.12  (-02) 
.37  (-03) 
.16  (-03) 
.80  (-04) 
.46  (-04) 


ejj  II  5 RATE 


.46 

(-02) 

1.95 

.20 

(-02) 

2.01 

.11 

(-02) 

2.00 

.72 

(-03) 

2.00 

.50 

(-03) 

2.00 

i 
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2 3 

Note  that  the  L -norm  of  the  error  for  the  spaces  S-(A  ) and 

0 B ,N 

A ^ 

S‘?(A  ) appears  to  be  smaller  than  that  for  the  more  heavily  graded 

U P , N 

3 4 

spaces  S_(Aq  „)  and  while  the  S-norm  of  the  error  is  larger. 

Mor'jover,  it  appears  that  the  rate  of  convergence  in  the  S-norm  has  been 
reduced  from  to  0(N~^^^).  A slight  modification  of  the  proof 

of  Theorem  6.1  shows  that,  in  general,  using  B > 6^  yields  this  rate  of 
convergence  in  the  S-norm. 

In  order  to  illustrate  the  result  of  Theorem  6.2,  we  computed  the 

2 

L -projections  of  u on  several  spline  subspaces  using  only  the  grading 

required  in  that  theorem,  6 = B^.  The  results  are  given  in  Tables 

2 

6. 7 - 6.9.  The  data  confirm  that,  with  B = B^,  the  L -projection  is 

k*'^-order  accurate  (in  the  L^-norm) ; it  also  seems  that  the  rate  of 

2 

convergence  in  the  S-norm  is  only  1/4  that  of  the  L -rate.  Moreover, 

2 

the  convergence  at  the  knots  is  only  at  1/2  the  L -rate  — the  opposite 
of  superconvergence! 
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N 

II  ejj  II  0 

RATE 

II  ejg  II  ^ 

RATE 

II  II  S 

RATE 

8 

.40  (-02) 

.82  (-02) 

.19  (+00) 

16 

.47  (-03) 

3.09 

.29  (-02) 

1.51 

.10  (+00) 

0.84 

24 

.14  (-03) 

3.02 

.16  (-02) 

1.50 

.76  (-01) 

0.77 

32 

.58  (-04) 

3.00 

.10  (-02) 

1.50 

.61  (-01) 

0.76 

48 

.17  (-04) 

2.99 

.55  (-03) 

1.50 

.45  (-01) 

0.75 

Table  6.7: 

2 

Error  in  L -projection 

on  8^(43 

n)  — 3 = 

N 

II  ®N  II  0 

RATE 

II  II  00^  A 

RATE 

II  e,  II  s 

RATE 

8 

.35  (-03) 

.27  (-02) 

.99  (-01) 

16 

.41  (-04) 

3.  10 

.67  (-03) 

2.04 

.49  (-01) 

1.00 

24 

.76  (-05) 

4. 13 

.29  (-03) 

2.00 

.33  (-01) 

1.00 

32 

.24  (-05) 

4.01 

.16  (-03) 

2.00 

.24  (-01) 

1.00 

48 

.50  (-06) 

3.84 

.74  (-04) 

2.00 

.17  (-01) 

1.00 

Table  6.8: 

2 

Error  in  L -projection 

^0(^4,n)  — 4 = B3 

N 

II  ejj  II  g 

RATE 

II  II  4 

RATE 

II  ej^  II  g 

RATE 

8 

.17  (-03) 

.85  (-03) 

.58  (-01) 

16 

.60  (-05) 

4.81 

.10  (-03) 

2.98 

.20  (-01) 

1.49 

24 

.67  (-06) 

5.42 

.32  (-04) 

2.99 

.11  (-01) 

1.49 

32 

.13  (-06) 

5.67 

.13  (-04) 

3.00 

.73  (-02) 

1.50 

48 

1 

.13  (-07) 

5.67 

.40  (-05) 

3.00 

.40  (-02) 

1.50 

2 

Table  6.9:  Error  in  L -projection  on  Sq(A^  — 6 = 


f 
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6.5  A Comment  on  the  Computational  Complexity  of  the  Algorithm 

Mesh  grading.  In  effect,  puts  more  of  the  computational  effort  Into 
the  region  where  u Is  badly  behaved.  The  price  one  pays  Is  the 
corresponding  decrease  In  effort  (and  accuracy)  elsewhere.  As  shown  In 
Section  4.3,  the  largest  Interval  In  a B-graded  mesh  Is  approximately 
B times  the  size  of  the  corresponding  Interval  In  a uniform  mesh. 

Thus,  to  obtain  the  same  accuracy  as  we  obtain  In  a nonsingular  problem 
using  a uniform  mesh,  we  might  have  to  take  a mesh  with  B times  as 
many  Intervals.  The  computational  complexity  of  the  numerical 
procedures  used  In  the  RRG  method  Is  0(N)  for  a mesh  with  N Intervals. 

3 

(To  set  up  the  linear  system  of  (2.10)  takes  0(k  ) work  In  each 
sub Interval.  The  order  of  the  system,  which  Is  equal  to  the  dimension 
of  the  spline  subspace,  grows  linearly  with  N;  since  It  Is  banded  with 
bandwidth  Independent  of  N,  Its  solution  requires  0(N)  work  and 
storage.)  Therefore,  this  method  requires  on  the  order  of  B times  as 
much  work  for  the  same  accuracy  as  In  the  nonsingular  case. 


1 


CHAPTER  7 
WEIGHTED  SPLINES 


7. 1 Introduction 

In  this  chapter  we  consider  approximations  of  the  form  x~°s(x), 
where  s(x)  is  a smooth  function  vanishing  at  0 and  1.  When  s(x)  is  a 
piecewise  polynomial,  we  call  these  functions  "weighted  splines".  We 
know,  from  Theorem  5.1,  that  u = x v,  where  v is  smooth.  Therefore,  we 
expect  that  x v will  be  a good  approximation  to  u when  v is  a good 
spline  approximation  to  v. 

In  Section  7.2,  we  define  weighted  splines  and  obtain  a useful 
bound  on  the  S-norm  of  a function  w,  vanishing  at  0,  in  terms  of  the 
L“-norm  of  D(x  w) . In  Section  7.3,  we  prove  error  bounds  for  a weighted 
spline  approximation  scheme  on  a weakly  6-graded  mesh  (1  ^ 3 < 2)  and 
use  this  result  to  obtain  optimal-order  error  bounds  for  the  RRG 
approximation.  As  the  approximation  v we  use  the  quasiinterpolant  of  v. 
In  contrast  to  the  8-grad Ing  used  in  Chapter  6,  here  6 remains  bounded 
above  for  fixed  k and  all  0 £ o < 1.  Moreover,  in  Chapter  6,  for  fixed 
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o,  8 -►  “>  as  k but  In  this  chapter,  for  fixed  o,  8 -►  1 as  k 

Numerical  results  are  included  in  Section  7.4, 


7,2  Weighted  Sp 1 ine  Subspaces 

Definition  7.1.  If  0 ^ c < 1 and  S^  is  a space  of  functions  defined  on 
1,  let 

S = {x“'’s  I s e S }. 
n , “ n 

G 

If  is  any  of  the  spaces  S^(4,^)  (respectively,  5^(4,^),  Sq(4)),  then 
we  call  the  elements  of  S^  ^ "weighted  splines"  and  denote  the  space 
^n  o Ws''^(o,4,^)  (respectively,  WSq(o,A,£),  WSq(o,4)). 

The  following  result  will  be  used  both  to  show  that  WSQ(a,A^^)  c g 
and  to  prove  error  bounds  in  Section  7.3. 


Lemma  7.1.  If  w(x)  = x v(x) , 0 < x £ 1,  where  v £ W^’  and  v(0)  = 0, 
then  w e S and 


w 


S(0,b) 


Dv 


L"(0,b)  all  0 < b < 1, 


with  Cq 


1+0 


Proof : By  L'Hopital's  rule. 


lim  w(x) 


lim 


v(x) 
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11. 


0_1 

since  x while  Dv(x)  is  bounded.  Thus,  w(0)  may  be  taken  to  be  0. 

Furthermore, 


S(0,b)  ■ ' 


x°  [x'^v  - Ox“°  -]^  dx 
0 * 


J"  x“°  [(Dv)^  - 2®(-)Dv  1-  dx. 

0 XX 


I “ 

Since  v(0)  = 0 and  v ^ W^’  , 


I - I = - 

X 

^ Dv(t)  dt 

I X I X 

0 

^ ll°''llL"(0,b) 


for  all  0 < b < I,  so  that 


Therefore 


X ''L"(0,b)  - " L“’(0,b) 


S(0,b)  ^ (I  + 2"  + J x-'’  dx 

= c^b^-'’  < « . 


- 87 


Corollary:  If  S_  c then  S _ c s. 

n 0 n,a 

Proof : If  V e W^’  and  w = x *^v,  then  w(l)  = v(l)  = 0.  The  preceedlng 
lemma  shows  that  w(0)  = 0 and  that  x°^^Dw  £ L^(I). 

G 


In  particular,  WSQ(a,A,£)  c g. 


7 . 3 We ighted  Spline  Approximation 

According  to  Theorem  5.1,  If  u e S'",  then  v = x^u  e We  now 

consider  using  u = x”°v  to  approximate  u,  where  v e F v is  the 

a 

quasiinterpolant  of  v« 


Theorem  7.1:  If  u e S'",  ra  > 0,  A = A with  B = 6.  e - — 

B,N  A 2(*— 1)  - o 

Ir 

is  any  incidence  vector,  then  there  exists  an  element  u £ WSq(o,A,^) 
satisfying 


, and  z 


l|u-u||g  < c^  ||d‘'"^(x'’du)  IIq, 

where  - mln(k,  m+2 ) and  c^  is  independent  of  u and  N. 

Proof : Let  v = F^v.  We  have  already  shown  that  if  = 0 and  » 

then  V e i.e.,  it  satisfies  the  boundary  conditions.  Moreover, 

by  Lemma  3.3,  if  we  also  take  = Q,  then 


D(v  - v) I . = 0. 
x=0 


This  result  enables  us  to  use  Lemma  2.1  to  bound  ||D(v  - v) 


L"(l^) 


which,  in  turn,  bounds  ||  u - u 


S(I^)’ 
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UNCLASSIFIED 
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We  first  consider  1^.  By  Theorem  5.1.  v « Thus. 

V - 9 € W^’  (Ip  and  v(0)  ■ 9(0)  ■ 0.  By  Lemma  7.1. 

" “"sap  ^ *=0  ‘'1'°  "l-(i,)* 

since  the  quasllnterpolant  Interpolates  Dv  at  0.  Lemma  2.1  applies,  and 

i "I'" 

Now  using  the  error  bound  (3.9)  for  the  quasllnterpolant  and  the  bound 
(4.9)(d), 


l|D(v  - 9)  < K h[/^  |ej‘-2 

< K C(6)‘-2  ||d‘v 

(By  assumption  1 = min(k,  nH-2)  ^2.)  Thus, 

"“■“"s(ip  - (co‘^C(B)‘"^)^  llD*v||p(ep 


(7.1) 


N 


-2(1-1) 


II  D^  II 


2 

l'(0^)* 


since  ,2(t-l)-o  .^-Bt2(l-1)-^)  .N-2(i-l). 

In  the  other  Intervals,  we  have,  by  (3.9)  and  the  bound  (4.9)(d), 

•1“  ■ “IILt  ^ ^ x"®[(D(v-9))^  - 2o(-5(=5)d(v-9)  + o2(^)2jdx 

i 


-a 

^ ’'i-i 


I|D(v-9)  llp(j  j +2o||D(v-9)  IIl2(j  j 


*■!_  1 II  " L*  (I  ) 


+ x^2^  ||v-9  II  2j 


L‘(ip 


(7.2) 
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- *1-1  ^ + o2(^)2j  1|d‘v1|J 


■ " 'x  ' " " L*  » 

i-l  *i-l  *'  '®i' 


But,  by  (4.9)(c), 


1 < M(Ab  .)  < C(B). 

'i-l  ‘'l-l 


and  by  (4. 9) (a). 


2(l-l)(a-l) 


-o  .2(1-1)  ^ B2(t-l)£^lJ_^ 

^ ~ (1-1)®° 


,-B  (2(Ul)-o] 


< 22(‘-1)  2®°  N-®  12(1-1)-°) 

< 4* 


(Here  we  have  used  the  identity 

2(B-l)(l-l)  - B(2(l-l)-o)  - 2(t-l)  + Bo  - Bo 


and  the  inequalities  Bo  £ B £ 2 and  1 2.)  Thus,  by  (7.2), 

Together  with  (7.1)  and  Theorem  5.1  (again  noting  that  at  most  2k  o£  the 
0^  overlap) , this  yields 

N 

llu  - u||2  £ I Hu  - u II  2 

1-1 

N 

£ max(c2,  Co)  i II  II  fS  x 

2 3 1-1  ^ <®i) 

£ 2k  max (c^,  I1D*^v|1q 

£ 2kr2  max(c2,  C3)  ||D^"^(x°Du)  ||  q. 


F 


2 

We  now  obtain  S-norm  and  L -norm  error  bounda  for  the  MIC 

approximation  to  u in  weighted  spline  subspaces.  The  mild  0-grading 

2 

required  in  Theorem  7.1  suffices  for  L -norm  error  bounds  of  optimal 
order. 


Theorem  7 . 2;  Let  u e be  the  RRG  approximation  to  the 

generalized  solution  u of  (4.1)  - (4.2),  where  6 ■■  B,  s . — and  z 

4 2(1-1 ) - 0 ~ 

is  any  incidence  vector.  If  p,  q,  and  f satisfy  the  hypotheses  of 

Theorem  5.2,  then 

Hu  - ullg  < jc^r  Ilf  IIj_2 

and 

llu  - uIIq  < (Ac^r)2  N'‘  Ilf  II  j_2, 

where  t = mln(k,  nri-2). 

Proof:  The  result  follows  from  Theorem  7.1,  the  bounds  on  the  solution 
in  terms  of  the  data  (Theorem  3.2),  and  Nitsche's  trick  (Theorem  2.6). 

□ 


I 

( 

! 

I 

i 

I > 


I 


7.4  Computational  Considerations  and  Numerical  Results 

As  a practical  matter,  the  method  of  this  chapter  seems  to  be 
superior  to  any  other  so  far  presented  for  solving  (4.1)  - (4.2), 
Including  the  techniques  of  Chapter  6.  The  rate  of  convergence  in  the 
theory  is  optimal.  No  restriction  need  be  made  on  the  smoothness  of  the 
piecewise  polynomials,  as  is  the  case  in  Chapter  8.  Furthermore,  the 
computationally  convenient  basis  of  B-splines,  used  for  computing  with 


4 
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the  spline  spaces  S may  be  modified  in  an  obvious  manner  to 

provide  an  equally  attractive  basis  for  US  (o,A,£).  Only  the  slightest 
8-gradlng  is  required,  and  our  numerical  results  show  that  the  error 
actually  obtained  is  markedly  less  for  the  weighted  spline  spaces  than 
for  those  of  Chapter  6.  The  reason  for  this  was  given  in  the  note  on 
computational  complexity  in  Section  6.5. 

Using  the  functions  - x ® basis,  we  are  faced  with  a 

serious  problem  when  forming  the  matrix  A by  numerical  quadrature.  The 
Integrands  in  the  inner  products  of  (2.10)  have  a singularity  at  the 
origin  of  the  form  x ° , Using,  say,  Gaussian  quadrature,  the  error  due 
to  this  singularity  is  so  great  (especially  in  the  first  Interval)  that 
the  rate  of  convergence  of  the  RRG  approximation  is  seriously  affected. 
Increasing  Che  number  of  quadrature  nodes  in  Che  interval  does  little 
good  (see  Table  7.1  for  an  example  using  weighted  cubic  splines  on  the 
problem  of  Section  6.4;  the  integrals  were  computed  using  a 5-point 
Gaussian  quadrature  in  each  interval) . 


8 

.56 

(-01) 

.13 

(+00) 

.25 

(+00) 

16 

.2A 

(-01) 

1.2A 

.58 

(-01) 

1. 16 

.13 

(+00) 

o.% 

i 

2A 

.19 

(-01) 

0.5A 

.A8 

(-01) 

0.A6 

.12 

(+00) 

0.25 

• 1 

32 

.16 

(-01) 

0.5A 

.A2 

(-01) 

0.A8 

.11 

(+00) 

0.25 

• 

AO 

. lA 

(-01) 

0.5A 

.38 

(-01) 

0.A9 

.10 

(+00) 

0.26 

Table  7.1:  Effect  of  Quadrature  Error 
Subspace;  WsJ(.5, i.09  - 


The  problem  is  solved  by  using  quadrature  formulae  which  are 
adapted  to  the  type  of  integrand  encountered  in  forming  A.  Thus,  we  use 
a formula  which  is  exact  for  Integrands  of  the  form 

i x"*"  p(x) 

i 

I for  p a polynomial  of  maximum  degree.  These  "weighted  Gaussian 

quadrature  formulae"  may  be  computed  using  an  algorithm  of  Golub  and 
Welsch  [Gl].  Finding  a k-point  formula  for  the  weight  function  w(x)  on 
the  Interval  (a,b)  requires  that  the  first  2k+l  moments  of  u, 
b 

/ xJ  (d(x)  dx,  0 ^ j £ 2k 
a 

be  computed.  The  procedure  then  takes  the  Cholesky  decomposition  of  a 
k+1  X k+1  matrix  and  solves  a k x k symmetric  tridiagonal  eigensystem. 
The  eigenvalues  and  the  first  components  of  the  eigenvectors  are  the 


nodes  and  weights  of  the  quadrature  rule. 


1C  uaed  in  every  interval  o£  the  nesh,  this  procedure  can  increase 
Che  cost  of  computing  k.  However,  it  suffices  to  chose  a point  x in 
(0,1)  (say  X * .05)  and  use  the  special  quadrature  rules  only  in 
intervals  which  intersect  (0,x).  Tables  7.2  - 7. A illustrate  Che  result 
of  this  procedure,  applied  to  the  problem  considered  in  Section  6. A. 

Computations  were  again  performed  in  double  precision  (5A  binary 
digits)  on  a PDP-10.  The  matrix  A and  vector  ^ were  computed  using 
k-1  - point  Gaussian  quadratures  in  each  interval  (weighted  or 
unweighted  depending  on  whether  or  not  Che  interval  intersects  (0,.0S)). 
The  errors  were  computed  using  k-M  - point  quadratures,  weighted  in  Che 
same  intervals.  Further  explanation  of  notaCional  details  is  given  in 


Section  6. A 


N " -,A 


.91  (-03) 
.13  (-03) 
.40  (-04) 
.17  (-04) 
.92  (-05) 


Table  7.2:  Error  In  WSq(.5,  Aj  — 1.14  - 


.17  (-01) 
.42  (-02) 
.18  (-02) 
.10  (-02) 
.67  (-03) 


10 

(-02) 

13 

(-03) 

37 

(-04) 

16 

(-04) 

80 

(-05) 

.42  (-05) 
.12  (-06) 
.16  (-07) 
.38  (-08) 
.12  (-08) 
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CHAPTER  8 


GENERALIZED  L-SPLINES 


8. 1 Introduction 

L-spllnes,  functions  which  are  elements  of  the  null  space  of  a 

differential  operator  L*L  In  each  Interval  of  a partition  A,  are  a 

natural  generalization  of  polynomial  splines  (see  (S5)).  L-spllne 

subspaces,  with  L » ''p(x)Du , were  used  In  the  RRG  method  for  the 

singular  problem  (4.1)  - (4.2)  by  Clarlet,  Natterer,  and  Varga  ICI),  vrtio 

obtained  L"-norm  error  bounds  which,  for  a uniform  mesh,  are  of  order 

h^~°.  They  noted  that  a suitably  chosen  non-uniform  mesh  might  Improve 

the  rate  of  convergence;  In  fact,  their  approximation  scheme  Is 

2 

second-order  accurate  If  a S-graded  mesh  Is  used  with  6 - (This  Is 

a special  case  of  the  main  result  of  this  chapter.) 

L-spllnes  were  subsequently  generalized  by  Crouzelx  and  Thomas 
[C4].  Given  a partition  A,  they  considered  functions  which.  In  the 
Interval  are  mapped  by  the  operator  L L Into  a predetermined 
f Inlte-dlmenslonal  space  of  functions  P^.  In  general,  this  space  can 
differ  from  one  Interval  to  another;  however,  for  the  purpose  of 
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obtaining  approximate  solutions  of  a differential  equation  in  which  the 

right-hand  side  is  smooth,  it  is  no  handicap  to  assume  that  is  the 

space  of  polynomials  of  some  given  degree.  Taking  L ■ /p(x)Ou  and 

the  space  of  polynomials  of  degree  < k-2,  Crouzeix  and  Thomas  showed 
2 

that  the  L -norm  of  the  error  in  the  RRC  approximation  to  the  solution 
of  (4.1)  - (4.2)  was  0(h'^~°) . They  considered  only  C®  subspaces, 
however,  and  their  results  do  not  apply  to  smoother  generalized 
L-splines. 


In  this  chapter,  we  redevelop  Crouzeix  and  Thomas's  theory. 

Because  we  consider  only  problems  of  the  form  (4.1)  - (4.2),  our 

presentation  is  more  elementary;  moreover,  we  have  simplified  the 

proofs  of  Lemmas  8.1  and  8.2.  We  then  show  that,  using  a g-graded  mesh 

with  B “ T-  ■»  the  error  is  of  order  N“*^, 

2~o 


8 . 2 Generalized  L-spllne  Approximation  in  s”* 

Let  A be  a given  partition  of  I,  and  let 

Sq(I^)  Ms  c S I s(Xi_^)  - s(x^)  - 0),  1 < i < N. 

Note  that  Sq(i^)  is  a closed  subspace  of  S.  The  following  result  is  an 
analogue  of  the  Rayleigh-Ritz  inequality  (cf.  Theorem  2.1)  and  is 
similar  to  Lemma  5 of  [C4]. 


Lemma 


8. 1:  For  all  v c ^ 

Xi  Xi 

/ v^  dx  £ C(lj^)^  / x°(Dv)2  dx. 


‘1-1 


‘i-1 


(8.1) 
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where 


(8.2) 


c(i.)2  - 


.2  x® 

l-l 


i > 1. 


4x2-°  i - 1. 


Proof ; For  i > 1,  by  the  Raylelgh-Ritz  inequality. 


f v2  dx  £ h2  / (Dv)^  dx 


< h2  x^j  f 


x°(Dv)2  dx. 


For  i ••  1,  an  integration  by  parts  shows  that 


1 A A 

f v2  dx  - xv2  - 2 / xv(x)Dv(x)  dx. 

0 'o  0 

Since  v(0)  - v(x^)  . q,  the  Integrated  term  vanishes.  Therefore,  by  the 
Cauchy-Schwarz  inequality, 

f v2  dx  £ 2 r / v2  dx'|^^2  (f  (xDv)2  dx'^^^^. 


Cancelling  the  common  factor  ||  v 
inequality  yields 


^ and  squaring  the  resulting 


*1  *1 
I dx  < 4 / x^(Dv)^  dx 


£ 4x2-°  ; x°(Dv)2  dx. 

^ 0 
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Let  k be  an  Integer,  k ^ 2.  A o-polynoilal  (of  degree  k-l ) la  a 
function  of  the  form 


Cq  + + 


+ c 


k-l-0 


k-l  = 


Note  that  s is  a o-polynomial  of  degree  k-l  If  and  only  If  0(x°D8)  is  a 
polynomial  of  degree  k-3  (the  0 polynomial  when  k - 2).  We  denote  by 
P(o,k)  the  space  of  o-polynomials  of  degree  k-l. 


Definition  8. 1;  A function  8(x)  c S is  said  to  be  a (o,k)-spline  with 
respect  to  A if  it  coincides  with  a o-polynomial  of  degree  k-l  in  each 
interval  of  A.  We  denote  the  set  of  all  (o.kl-splines  with  respect  to  A 
by  S^(A). 

Q 


Clearly,  If  s e S^(4),  then  there  exist  polynomials  pj,  Pf/  of 

degree  k-3  such  that  Dlx^Ds)  - p^  in  Ij^,  1 £ 1 < N.  Integrating  by 
parts,  we  have  that 


(8.3) 


Xi  Xi 

/ x®DsDv  dx  - / -0(x‘’Ds)v  dx 

’‘i-1  ’‘i-1 

’‘i 

“ “Pi^  dx,  for  all  v * Sgd^) , 1 £ 1 1.  N. 


The  converse  also  holds.  Let  denote  the  restriction  of  s to  1^^. 

If  s e S*^  satisfies  (8.3),  then  D(x“ds^^P  . (i.e.,  s^^j  is  a 

o-polynomial)  in  each  interval  g^nce  Sgdi)  is  dense  in  L^d^). 


I 
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It  happens  chat  Che  orthogonal  projection  of  a function  u on 

a 

with  respect  to  ^ Interpolates  u at  the  luiots.  This  result  Is 

due  CO  Crouzeix  and  Thomas  1C4]. 


Lemma  8.2:  Given  u £ S,  there  exists  u £ S^(a)  satisfying 


(8.4) 


u(x^)  - u(x^) , 0 < I < N 


(8.5)  II  u - u 11  - inf  II  u - s 11  . 1 < 1 < N. 

^'^i^  8£P(o.W)  ' i'  “ “ 

Proof:  Suppose  s £ : {s  « S | D(x°Ds^j^j)  £ L^(l£)}  is  orthogonal 

to  the  subspace  Sq(i^)  with  respect  to  the  inner  product  (•♦•)s(i  j* 

i . e . . 


f x"DvDs  dx  » 0 for  all  v £ S.(l  ) 


[early,  s satisfies  (8.3)  with  p^  » 0.  Therefore,  is 


a-polynomlal : s^^^  £ P(o,2)  = P(o,k.) 


We  denote  by  X"*"  the  space  of  functions  orthogonal  (with  respect  to 
(.,.)s(i  j)  to  the  set  X.  We  have  just  shown  that 

s £ (S°(I^)  n Sgdi)^)  ®(i)  e P(<J,2)  «=P(o,k), 


(S°(I^)  n Sq(I^)^)  = P(a,k). 

Thus  P(o,k)‘*'  c (S®(I^)  n Sq(I^) Moreover,  since  S^(l^)  is  dense  in 

S(Ij^),  (s°(Ii)  n Sodi)"^)  ^Odi)"^*  ®"y*^hing 

orthogonal  to  (S®(1^)  n Sudj^)"*")  is  orthogonal  to  Sq(Ij^)'*'.  Hence, 
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P(o.lc)^  c Sq(1^)^^  - Sq(i^), 
since  Sq(I^)  is  a closed  subpsace  (see  Schecter  ISIJ). 

Now  let  t be  the  orthogonal  projection  of  u on  P(a,k).  Since 
u - t e P(o,k)'*’,  it  follows  that  u - t c ^^(1^^);  l.e.,  t interpolates 
u.  It  is  now  clear  how  to  construct  u.  We  simply  piece  together  the 
best  o-polynomial  approximations  in  each  interval;  the  resulting 
function  clearly  satisfies  (8. A)  and  (8.5),  is  continuous  and  satisfies 
the  boundary  conditions. 

□ 

We  now  obtain  error  bounds  for  u - u.  This  result  and  its 
corollary  are  due  to  Crouzeix  and  Thomas  [CA]. 

0 ..  k, 

Lemma  8.3:  Let  u e S and  u e be  the  approximation  defined  in 

Lemma  8.2.  Then  for  all  1 < i £ N, 

(8.6)  ||u-u||  . < C(I.)  inf  1|  D(x‘’du)  - p ||  , j .. 

^^^i^  - P€Pj^_2  ' 

(8.7)  ||u-u||  j < (C(I,))^  inf  ||  D(x“Du)  - p . 

^ i^  P^\-2 

Proof:  Let  p e ^^-2  arbitrary.  Let  Vp  ^ P(n,k.)  be  the  unique 
solution  to  the  boundary-value  problem 

D(x‘’dv)  “ p,  < X < x^, 

v(x^_^)  - u(x^_j),  v(x^)  - u(x^). 
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Obviously  u - Vp  c SQ(Ij^).  Integrating  by  parts,  using  the  boundary 
conditions  on  the  Cauchy-Schwarz  Inequality,  and  (8.1), 


P "s(ip 


f -(u  - V )[d(x°Du)  - pj  dx 


< U - V 


< C(I^)  Ilu-Vplis^^j  I1D(x"du)  - 


“■''plls(l^)  < Cdi)  IlDlx^Du)  - p 


But  u Is  the  best  approximation  to  u In  the  norm  ||  . 


S(i  ),  so 


“-""sdi,  - ■''p'lscij, 


< C(I  ) Inf  ||D(x'’Du)  -pH,,.,,, 


which  proves  Inequality  (8.6).  Since  u - u e Sq(i^),  inequality  (8.7) 


follows  from  (8.1)  and  (8.6). 


0 ~ k. 

Corollary;  Let  u e S and  u c be  the  approximation  defined  In 

Lemma  8.2.  Then 


(8.8) 


u - u |l  £ max  {C(I  ) } r(u) 
KKN 


(8.9) 


u - u II  < max  {(C(I  ))'‘}  R(u), 
KKN  ^ 
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where 


c 


(8.10)  R(u)  = 


N 

I 


inf  ||D(x°Du) 

V‘'‘  ''Vz 


1/2 


We  now  bound  the  quantities  C(l^)  and  R(u)  in  terms  of  N 


Sobolev  norms  of  D(x  Du) . 


Lemma  8.4;  If  A “ A „ with  B “ — , then  the  constants  C(I. 

PfW  z— 0 1 

in  Lemma  8.1  satisfy 


(8.11) 


,-l 


C(I^)  < 2N  \ 1 < i < N. 


Proof;  By  (8.2)  (for  i = 1), 


C(l^)2  - 4x^"° 


4N-^<2-a)  . ^^-2^ 


By  (8.2)  (for  i > 1 ) and  (4. 9) (a), 

, 2 


C(I^)‘ 


2 


IT  X 


i-1 
2 ,2(B-1) 


e i 


-B(2-o) 


- (1-1 


Since  6(2-0)  - 2,  2(6-1)  - 6o,  and  6a  < 6 < 2, 


b2  ^2(b-1) 


(i-1) 


6o 


< < 16, 


and 


cdi)^  1 "4  ^ ^ i 


□ 


and  Che 


) defined 
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Lemma  8._5  (Clarlet  and  Raviart  [02]):  If  u c s”*,  m ^ 0,  then 

(8.12)  R(u)  < h‘‘2  IID^-^Cx^Du)  IIq, 

where  R(u)  is  given  by  (8.10)  and  i = min(nH-2,  k) . 

Corollary:  Let  u ^ S®  and  u ^ Sg(A^  ^)  be  the  approximation  definea  in 

2 

Lemma  8.2  with  6 » '2~o' 

(8.13)  llu-u||g  < 26''“2  11  d‘'"\x‘’du)  11  q 

(8.14)  11u-u11q  £ 48*’“^  N *■  11  D*'  ^ (x'^Du)  11  Q. 

Proof ; Use  (8.11)  and  (8.12)  to  bound  the  right-hand  sides  of  (8.8)  and 
(8.9). 

□ 


Crouzeix  and  Thomas  obtain  error  bounds  for  the  RRG  approximation  u 

If 

to  u in  S^(A)  by  showing  that  u - u is  small.  We  prefer,  however,  to 
use  Nitsche's  trick  (Theorem  2*6). 


Theorem  8.1:  Let  u ^ S^(Aq  ) be  the  RRG  approximation  to  the 
o p 

2 

generalized  solution  u of  (4.1)  - (4.2),  where  ® “ 2-o‘  If  p,  q,  and  f 
satisfy  the  hypotheses  of  Theorem  5.2,  then 


l|u-ullg  < (2  ^ llflU^2 

Hu  - ullg  < (2A  B^-2r)2  Ilf  11i_2. 


where  - min(nri-2,  k) . 


PART  III 


CHAPTER  9 

SPHERICALLY  SYMMETRIC  PROBLEMS 


9. 1 Introduction 

In  Part  III,  we  consider  the  numerical  solution  of  spherically 
symmetric,  elliptic  partial  differential  equations.  There  are  numerous 
applications  in  engineering  and  the  sciences  in  which  the  solution  of  a 
spherically  symmetric,  elliptic  equation  is  desired  (see  the  references 
in  [R3]).  Since  all  the  functions  involved  are  spherically  symmetric 
(that  is,  they  depend  only  on  distance  from  the  center  of  the  domain), 
the  problem  can  be  replaced  by  an  equivalent  two-point  boundary  value 
problem. 

When  an  n-dimensional  problem  is  so  simplified,  the  resulting 
problem  has  the  form 

- D(r””V(r)Du)  + r""^(r)u  = r'^'^fCr),  0 < r < 1, 

u(l)  = 0,  Du(0)  “ 0 (alternatively,  u(0)  finite). 


In  marked  contrast  to  the  singular  problems  of  Part  II,  the  solution  is 
smooth  despite  the  singularity  in  the  equation.  It  should  therefore  be 
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possible  to  approximate  the  solution  accurately  using  the 
Raylelgh-Rltz-Galerkln  method  with  a piecewise  polynomial  subspace, 
without  any  sort  of  mesh  grading.  In  Chapter  10,  we  obtain 
optimal-order  error  bounds,  showing  that  this  procedure  is  theoretically 
well-founded.  Instead  of  the  usual  one-dlmenslonal  Sobolev  norms,  we 
use  norms  which  are  appropriate  to  the  original  n-dlmenslonal  setting  of 
the  problem. 

In  2 and  3 dimensions,  Russell  and  Shamplne  [R3j  proved  error 
bounds  for  approximation  procedures  specially  designed  to  deal  with  the 
apparent  singularity  at  the  origin.  In  particular,  they  treated 
collocation  in  which  the  basis  is  augmented  by  singular  basis  functions, 
RRG  using  singular  patch  bases  (L-splines),  and  a finite-difference 
scheme  of  Jamet  [Ji]  designed  to  handle  the  singularity.  Crouzeix  and 
Thomas  IC4]  and  Reddlen  [Rl]  considered  this  problem  as  part  of  a wider 
class  and  obtained  similar  results  for  subspaces  which  Include  singular 
basis  functions. 

Dupont  and  Wahlbin  [D6]  and  Jesperson  [J3]  analyzed  the  RRG 
procedure  with  piecewise  polynomials,  and  obtained  error  bounds  of 
optimal  order  in  the  one-dimensional  Sobolev  norms.  The  import  of  their 
results,  together  with  those  of  these  two  chapters,  is  that  no  special 
measures  are  required  for  this  problem:  the  Rayleigh-Ritz-Galerkin 
method  using  high-order  piecewise  polynomial  spaces  on  a uniform  mesh  is 
a highly  effective  numerical  method. 
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Section  9.2  Is  a summary  of  the  variational  form  of  the  problem  and 
the  properties  of  its  solutions.  Bounds  on  solutions  of  Hermite 
Interpolation  problems,  required  for  the  error  bounds  of  Chapter  10,  are 
proved  In  Section  9.3. 


9. 2 Spherically  Symmetric  Elliptic  Problems 

Let  B(b)  be  the  open  ball  of  radius  b in  r”,  i.e., 

B(b)  = {x  e r"  I r(x)  < b}, 

2 2 2 

_ Xj^  + X2  + . . . + , and  let  B = B(l)  be  the  open  unit 

ball.  We  say  that  a function  V(x)  defined  on  B is  spherically  symmetric 
if  it  depends  only  on  distance  r(x)  from  the  origin.  If  V is 
spherically  symmetric,  we  call  the  function  v(r)  such  that 
v(r(jc))  = V (x)  the  radial  part  of  V. 


Consider  the  spherically  symmetric,  elliptic  partial  differential 
equation 

(9.1)  - V(p(r(x))VU)  + q(r(x))U  = f(r(x)),  x c B, 

(9.2)  U(x)  - 0,  X e 3B, 


where  p,  q,  and  f are  given  functions  defined  on  [0,1].  We  assume  that 


(9.3) 


a) 

p(r(x) ) 

e C^B), 

b) 

p(r(x)) 

- '’min  > ° all  x ^ 

c) 

q(r(x)) 

« C(B), 

d) 

q(r(x)) 

^0  forallx^B* 
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Since  the  data  of  the  problem  are  spherically  symmetric,  an  obvious 
symmetry  argument  shows  that  U is  too.  (A  change  in  coordinate  systems 
by  rotation  around  any  axis  passing  through  the  origin  leaves  the 
problem,  and  hence  its  solution,  unchanged.) 

By  transforming  equation  (9.1)  to  spherical  coordinates  and  setting 
to  0 all  partial  derivatives  with  respect  to  angle,  one  obtains  a 
singular  two-point  boundary  value  problem  for  the  radial  part  u of  U: 

(9. A)  - D(r"  ^p(r)Du)  + r”  ^q(r)u  - r”  ^f{r),  0 < r < 1, 

(9.5)  u(l)  “ 0,  Du(0)  “ 0 (alternatively,  u(0)  finite). 

When  n - 3,  the  change  of  variables  v • xu  results  in  the  nonsingular 
problem 

- D^v(r)  + q(r)v(r)  « rf(r),  0 < r < 1, 
v(0)  =■  v(l)  = 0. 

Unfortunately,  we  know  of  no  such  trick  when  n V 3! 

Another  one-dimensional  analogue  of  (9.1)  - (9.2)  is 
(9. A')  - D(|xr“^p(|x|)Du)  + |x|‘‘”^q(|x|)u  - |x|  ""^f  ( |x|  ) , -1  < x < 1, 

(9.5')  u(-l)  - u(l)  “ 0. 

The  relation  to  (9. A)  - (9.5)  is  that  the  solution  u(r)  of  the  former 

problem  is  the  restriction  to  [0,1]  of  the  solution  of  the  latter. 

i 

t 


f. 


- 108  - 


Suppose  that  G Is  spherically  symmetric  and  that  g Is  the  radial 
part  of  G.  Using  an  (r,j)  coordinate  system  In  which  r - r(x)  Is  given 
by  (9.1)  and  ^ represents  an  n-1  dimensional  vector  of  angles,  we  have 

b 

G dx  - / r"  G(r,0^)  d9.dr 

B(b)  0 

b , 

- n ; r”  g(r)  dr, 

" 0 

where  B Is  the  area  of  the  unit  hypersphere  In  IR*'  IM3] . This 
n 

observation  motivates  the  following  definition. 


Definition  9.1;  For  real-valued  functions  f,g  defined  on  (0,b),  let 


n-1. 


(f,g)„,^,  = / r“  "f(r)g(r)  dr 


B(b) 


and 


B(b) 


= (f 


□ 


Definition  9.2:  Let  J™(b)  (respectively,  JQ(b))  denote  the  closure  of 


the  C“  functions  with  all  odd  derivatives  vanishing  at  0 (respectively, 
the  C"  functions  which  vanish  in  a neighborhood  of  b with  all  odd 
derivatives  vanishing  at  0)  with  respect  to  the  norm 


m,B(b) 


m 

Z 


D^f 


1/2 


J-0 


B(b) 


□ 


Note  that  J°'(b)  (respectively,  J™(b))  may  be  identified  with  the 


J 
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restriction  of  the  spherically  symmetric  functions  in  H'''(B(b)) 
(respectively,  H“(B(b))  to  any  line  segment  from  the  origin  to  the 
boundary  of  B(b) . 

Our  assumptions  on  p and  q Imply  that  there  exist  positive 
constants  X and  A such  that 

(9.6)  < •<“.'‘>8(6)  i Mb) 

for  all  u € J^(b),  where 

a(u,v)  = / (pDuDv  + quv]  dr. 

B(b)  Q 

The  results  of  Sections  2.2  and  2.3  apply  to  (9.4)  - (9.5).  J®  is  the 
underlying  Hilbert  space,  and  the  space  S of  admissible  functions, 
with  norm  ||  v 1|  ^ - ||  Dv  1|  g. 

Our  approximation-theoretic  results  in  the  spaces  will  rely  on 
the  following  basic  result. 

Lemma  9. 1 (Friedrich's  inequality  IM3]):  If  v e jQ(b) , then 

(9.7)  llvlljj^,  < b IIDvIl^,^,. 

Inequalities  (9.6)  and  (9.7)  show  that  the  bilinear  form  a(  , ) is 
positive  definite  over  the  space  J^.  Thus,  for  each  f t there 

exists  a unique  u ^ J^,  the  generalized  solution  of  (9.4)  - (9.5),  such 
that 

a(u,v)g  . (f,v)g 


for  all  V c Jq. 
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Regularity  results  for  the  problem  (9.1)  - (9.2)  are  well-known, 
and  corresponding  results  for  (9.4)  - (9.5)  follow  immediately. 

Lemma  9 . 2 [F3);  There  exists  a constant  r such  that  for  all  m > 2,  if 
f € then  the  generalized  solution  u c J®  and 


(9.8) 


R i r II  9 R • 

m,B  ~ ra-2,B 


We  now  state  a variant  of  the  Sobolev  lemma.  Let 

”0  ’ = [ij 

Lemma  9 . 3 [F3]:  There  exists  a positive  constant  such  that,  if 
u e j“(b) , m > m then 


(9.9)  (u(0))2  < c2  I IIdJuII* 


''  j=0 


B(b) 


When  bounds  on  the  value  at  a point  (other  than  0)  of  a function  in 


J™(b)  are  desired,  the  Sobolev  lemma  is  not  the  best  result.  The 
following  stronger  result  shows  that  j“(b)  behaves  more  like  a 
one-dimensional  Sobolev  space.  Let  in  be  the  largest  odd  integer  less 
than  or  equal  to  m - m^,  i.e.. 


m - m^  - 1 (m  - m^  even) 


m 


"0 

m - niQ 


(m  - niQ  odd) 
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Lemma  9. A;  If  v « J (b)  then,  for  all  0 < a < b. 


a)  V £ H“(a,b); 


b)  V 


£ C“"^[a.b]; 


c)  11 


^(a.b) 

< cj  a^~"  l(b-a)"MlD^v  II  + (b-a)  ||  ||  ] 


0 ^ j £ m-1 ; 


d)  V has  m - m^  continuous  derivatives  at  0,  and  its  odd 
derivatives  vanish  there,  i.e.. 


(9.10) 


D v(0)  - 0,  i “ 1,  3,  ...,  m; 


e)  if  in  addition,  v e J™,  then 


D v(b)  - 0,  0 < i < m-1, 


Proof:  First,  note  that 

2 


m b 


llg  II 


n-1  ,„i  ,2 


m,B(b) 


^ Z / X (D  g)  dx 
i-0  a 


, m b - 

> a"“^  Z S (D^g)"^  dx 

i“0  a 


n-1  11  I,  2 

-a  II  8 II  „ 

H (a,b) 

Thus,  if  a sequence  of  functions  {v  } converges  to  v with  respect  to  the 

n 

J™(b)-nonn,  it  also  converges  in  the  H”(a,b)-norm.  This  immediately 


Implies  (a).  Next,  (b)  follows  from  (a)  by  the  one-dimensional  Sobolev 
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inequality  (Theorem  2.3).  Moreover, 


" " L*(a,b) 

< cj  [(b-a)-MlDJv||2,^^^^^  + (b-a) 


< cj  a^""  [(b-a)"^  IlD^v  + (b-a)  |1  0^"^^  ||  J 


which  proves  (c). 


By  the  inequality  (9.9),  convergence  in  the  j'”(b)-norm  implies 
convergence  at  zero  of  derivatives  up  to  order  m - m^.  Since  v is  the 
limit  in  this  norm  of  a sequence  of  smooth  function  satisfying  (9.10), 
v does  too,  proving  (d).  Finally,  (e)  follows  by  the  same  argument, 
using  part  (c)  rather  than  Lemma  9.3. 


□ 


9. 3 On  Hermite  Interpolation 

In  Chapter  10,  we  will  obtain  an  error  bound  involving  the 
(spherical)  norms  of  a polynomial  characterized  by  interpolation 
conditions  on  the  polynomial  and  its  derivatives  at  two  points.  In  this 
section  we  consider  bounds  on  the  solutions  of  such  problems. 

Let  k be  a positive  integer,  e^,  ej^  be  integers  such  that 
e^  + = k,  and  b > 0 be  given.  The  two-point  Hermite  interpolation 

problem  is  to  determine  a polynomial  Pj^(x)  of  degree  k-1  such  that 


hi. 
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(9.11) 


o^Pklo)  - 0<J<.„, 

I>’p^(b)  - y5^\  0 < j < ej. 


It  is  known  that  (9.11)  has  a unique  solution  for  arbitrary  data  y 
(S2]. 


(j) 


We  shall  need  bounds  on  the  norm  of  the  solution  of  an  Hennite 
Interpolation  problem  In  terms  of  the  size  of  the  data. 


Lemma  9.5;  There  exists  a constant  - C(eQ,ej)  Independent  of  b and 
y^J^  such  that 


(9.12) 


"“^Pk"B(b)  - ‘^e  ~ t Z bJ  \y[^^\ 

^ ® 1-0  j-0  ^ 


for  all  0 ^ I £ k,  where  p^  Is  the  polynomial  which  satisfies  (9.11) 
with  data  y|^\ 


Proof ; Let  z^^^  = b"^ 
that 


y^  and  be  the  polynomial  of  degree  k-1  such 


(9.13) 


Clearly,  Pj^(x) 


D^qj^(O)  = Zq^^  , 0 < j < eQ, 

D^qj^d)  - zJJ\  0 < j < e^. 


q|^('g’) . The  coefficients  £ 


Mx  - z. 


(Yj^)  of  qj^  satisfy 


where  M = M(e^,e2)  is  a nonsingular  k x k matrix  and  £ Is  the  data  of 
(9.13),  in  some  given  order  [S2].  Hence, 
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111  111  < IIm'MIj  IIzIIj, 


where  H^H  = I |z  | is  the  usual  vector  1-nonn  and  1|  ||  the 

i ‘ 

corresponding  matrix  norm.  Therefore 


D*'p 


k " B(b) 


"B(b) 


'[^i  "B(b) 


< V llD^(^)SllB(b) 


/ 


< II L II  . max 
i<i<k 


; x""^  [*'n  (i-j)  b"^  x^"*’]^  dx 

0 j=0 


C b^2  - lll_||  ^ 


< C 11  M~^  11  ^ b^2  " II  z 


r 


CHAPTER  10 


SPHERICAL  SPLINE  APPROXIMATION  THEORY 


10.1  Introduction 


If  V c H (B)  is  a spherically  symmetric  function,  we  propose  to 


approximate  it  with  smooth  spherically  symmetric  functions.  Let  v c J 


be  the  radial  part  of  V and  v 6 J be  a spline  approximation  to  v;  we 


take  as  our  approximation  to  V the  function  V whose  radial  part  is  9. 


As  a measure  of  the  error  (V  -V),  we  use  the  norm  ||dJ(v  - 9) 


B(b)‘ 


Spline  subpaces  SS^(A,^)  of  J™  and  SS^(A,^)  of  J™  are  defined 


in  Section  10.2.  First,  we  construct  spline  spaces  ES  (a',^')  and 


ESq(A',^')  whose  elements  are  even  functions  on  (-1,1),  requiring 


sufficient  smoothness  at  0 for  odd  derivatives  of  the  appropriate  orders 


V 

to  vanish.  The  restriction  to  [0,1]  of  s e ES  (A',z')  is  an  element  of 


J . we  then  define  the  space  of  spherically  symmetric  splines. 


SS  (A,£),  to  be  the  space  of  all  such  functions. 


Section  10.3  deals  with  approximation  by  spherical  splines.  As  the 


approximation  v,  we  use  a quasiinterpolant  In  Theorem  10.1,  we 


show  that  for  v e J , 
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I|dJ(v  - F^v)  llg  < C h“"J  I1d“v  II  g 

for  a partition  with  mesh-length  h,  where  the  constant  C is  independent 
of  V and  h.  We  obtain  similar  error  bounds  for  the  Rayleigh-Rltz- 
Galerkin  approximation  in  Theorem  10.2.  Numerical  results  are  included 
in  Section  10.4. 


10.2  Spline  Sub  spaces  of  J°* 

k.  ni 

We  shall  construct  a space  of  splines  SS  (A,^)  c j (respectively, 
SSq(A,^)  J™  n Jq).  In  view  of  Lemma  9.4,  we  must  force  the 
derivatives  of  orders  1,  3,  ...,  m of  such  splines  to  vanish  at  0.  Let 
k be  a positive  Integer,  A a partition  of  [0,1],  and  z an  incidence 
vector.  Define  the  symmetric  partition  a'  by 


A':-l  - x_^  < < ...  < Xq  - 0 < xj  < ...  < Xj^  = 1, 

where  ^ ^ 1 £ i £ N;  and  the  symmetric  incidence  vector  z'  by 


~ ° ^^-N+r  ^0’  •••’  ^N-P' 


where  ^ ^ = z^,  1 £ i £ N-1  (Zq  is  as  yet  arbitrary).  The  resulting 
knot  vector  £(A',£')  is  also  symmetric. 


Definition  10.1:  Let  ES^(A',^')  (respectively,  ESq(A',^'))  be  the  space 
of  even  spline  functions  of  order  k with  respect  to  A'  and 
(respectively,  the  even  spline  functions  of  order  k with  respect  to  A' 
and  z'  which  vanish  at  -1  and  1).  Let  SS*^(A,^)  (respectively,  SSq(A,2)) 
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be  the  space  of  functions  on  [0,1]  which  coincide  with  the  restriction 

k W. 

to  [0,1]  of  an  element  of  ES  ( A' ) (respectively,  ESq(A',£')). 


To  motivate  what  follows,  we  shall  outline  the  proof  of  error 
bounds  given  in  Section  10.3.  Let  be  quasiinterpolation  points  for 

which,  like  a'  and  , are  symmetric  about  0;  by  the  symmetry 
of  _t(A^,£^)  and  the  quasiinterpolant  maps  even  functions  into  even 

functions,  i.e.,  into  ES  (a',^').  We  shall  also  require  that  the  be 
bounded  away  from  0,  because  we  want  to  use  Lemma  9.4  (part  (c))  to 
bound  the  value  of  a function  at  tj . To  enforce  both  this  requirement 

V 

and  symmetry,  must  be  chosen  so  that  the  dimension  of  S (a'.^') 

(which  is  the  number  of  quasiinterpolation  points)  is  even. 


Another  constraint  on  is  that  we  shall  not  be  able  to  prove 
error  bounds  for  even  splines  which  are  forced  to  be  too  smooth  at  0. 
In  particular,  we  want  z^  to  be  sufficiently  large  so  that  elements  of 
ES  (a',£')  aren't  required  to  have  m + 1 continuous  derivatives  at  0. 


max(k-l-m,  1)  (k  odd) 


max(k-l-m,  2)  (k  even) 


il)  the  dimension  of  S (A',^')  is  an  even  number. 
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ill)  elements  of  S^(A',^')  and  ES^(A',^')  are  not  required 
to  have  a continuous  derivative  of  order  m + 1. 


Indeed,  (i)  is  obvious.  For  (ii),  note  that 


k 

d = dim(s‘^(A' ,z'))  - k,+  r z 

i=l-N 


= k + z + 2 2 z ; 

i=l 

d is  even  since  k + z^  is  even.  Finally,  the  elements  of  S and 

ES^(A',z')  have 


k - 1 - z , 


min(m,  k-3)  (k  even) 


min(m,  k-2)  (k  odd) 


continuous  derivatives  at  0,  which  proves  (iii). 


Because  the  par tition  _t (a' ,£' ) is  symmetric,  the  B-spline  basis 
functions  have  the  symmetry  property 


(10. 1)  N (-x)  “ N , , (x) , -l<x<l,  l£i<d. 


The  functions 


®j  " ^j,k‘^Vi-j,k 


. 1 1 J 1 7 . 


are  a basis  for  ES  (A',^'),  and  their  restrictions  to  [0,1]  are  a basis 

k.  k 

for  SS  (A,£).  (Bj^  must  be  deleted  from  bases  for  ESq(A',£')  and 

SSq(A,£).)  Clearly,  SSq(A,£)  is  a subspace  of  j”  n 
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10.3  Error  Bounds  for  the  Quasllnterpolant 

We  first  consider  approximation  by  polynomials  in  a neighborhood  of 
the  origin. 


Lemma  10.1:  If  v e J (b),  m ^ 1,  then  there  exists  a nolvnomlal  T.  v of 

b 

degree  m-1  satisfying 

(10.2)  llB(b)  - "““''I'BCb)’ 


Proof : According  to  Lemma  9.4,  v has  m-1  continuous  derivatives  at  b. 
Let  Tj^v  be  the  first  m terms  of  the  Taylor  series  for  v at  b,  i.e., 


T v(x)  = v(b)  + Dv(b) (x-b)  + . 

D 


D™  ^v(b) 


f^(x  - b) 


m-1 


We  prove  (10.2)  inductively,  starting  with  the  case  j * m and  proceeding 
downwards.  For  j = m,  (10.2)  holds  trivially.  By  the  interpolation 
conditions  which  define  T v,  D-^lv  - T v)  t ji(b)  for  0 < j < m. 
Therefore,  by  (9.7)  and  the  inductive  hypothesis. 


< b >l°"''llB(b)^- 

0 
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Lemma  10.2:  If  v € J (b)  , m ^ 0,  then  there  exists  a jxjlynomlal  of 
degree  m-1  such  that 

(10.3)  - S^v)  ll‘^“''llB(b) 

(where  is  independent  of  v and  b) , and  the  odd  derivatives  up  to 

order  m of  S,  v vanish  at  zero,  i.e., 
b 

(10.4)  D^S,v(0)  =0,  i = 1,  3,  ...,  m. 

b 

Proof:  If  m < m^,  then  (10.4)  is  vacuous:  choosing  S,  v = T,  v gives  the 
U b b 

desired  result.  Now  assume  that  m ^ m^.  According  to  Lemma  9.4,  v has 

m-1  continuous  derivatives  at  b and  ni“™Q  continuous  derivatives  at  0. 

Let  S V be  defined  by  the  m-1  conditions 
b 0 

D^Sj^v(b)  = D\(b),  i = 0,  1,  ...,  m^  - 2, 
and  the  m - m^  + 1 conditions 

(10.5)  D^S^v(O)  “ D^v(O),  i = 0,  1,  2,  ...,  m - m^. 

Sj^v  is  well-defined  since  it  is  the  solution  of  an  Hermite  interpolation 
problem.  Moreover,  (10.4)  is  satisfied,  because  by  Lemma  9.4  part  (d) , 
the  odd  derivatives  up  to  order  m of  v vanish  at  0. 

We  shall  bound  HD-'lv  - S.v)  lln...  by  showing  that 

b Dv  b) 

IlD^dj^v  - Sj^v)  II  is  of  the  same  size  as  ||D^(v  - T^^v)  llg^j,)*  But 

if  E V = Tv  - S V,  then 
b b b 

(10.6)  D^E,^v(b)  - 0,  i - 0,  1,  ...,  - 2, 

b U 


121 


L 


and  by  the  conditions  (10.5), 


(10.7)  D^E^v(O)  = - v)(0),  1-0,  1,  2 m - ihq. 


Since  solves  the  Hermlte  Interpolation  problem  (10.6)  - (10.7), 


bounds  on  E v follow  from  Lemma  9.5.  We  have 
b 


m-m„ 


IID^Ej^v  II  B(b)  < Z |D^(v  - TbV)(0)|, 


1-0 


and  by  Sobolev's  lemma  and  (10.2), 

m-iDo 


II  D^E^v  II 


o/■u^  < C Z b 

® 1=0 


(n/2)-j+l 


/ 


m^ 


"n  II^^b) 


1/2 


m- 


”o 


< C Z b 
® 1=0 


(n/2)-j+l 


II  d“v  II 


V 


B(b)J 


CeCn(nJ-mo+l  )/(Oo+l ) b“  ^ H II  B(  b)  ‘ 


□ 


We  require  the  following  bounds  on  weighted  norms  of  derivatives  of 
the  B-spllnes. 


Lemma  10. 3:  For  1 > 0 and  each  Integer  £ ^ 0, 


||^(n-l)/2  ^ 11^^ 


n/2 


< c. 


j,k"L^(I^)  - '-2  ^£  ’ 

1 


where  c^  depends  only  on  n and  £. 


Proof : By  Lemma  3.1, 
ii.(n-l)/2 


D N,  II 


j,k  "LdlJ 


; X . (x))^  dx 

X ^ 

1-1 
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I 2 * 1 

< IId‘n  ir 


< (B  h — 
— t 1 n 


□ 


We  now  show  that  to  every  v e J™  n there  corresponds  an 

th  ^ k. 

m -order  accurate  approximation  v in  SSq(a,^).  We  first  extend  v to 

X e (-1,1)  by  setting 

v(-x)  = v(x) , 0 < X < 1 . 

k 

Next,  we  show  that  the  quasiinterpolant  c ESr,(A',z')  is  m -order 

A u — 

accurate;  its  restriction  to  [0,1]  will  be  the  desired  approximation  in 
SsJ(A,z). 


Theorem  10.1:  There  exists  a positive  constant  c^  which  depends  on  the 
local  mesh  ratio  M(a)  such  that,  if  v e n J^,  m £ k,  then  there 


k 

exists  V e SSq(a,£)  such  that 

,1 


D'Cv  - ^)  II  g < h™  *■  11  D^v  II  g 


for  all  0 < I < m. 


Proof : Let  Be  quasiinterpolation  points  for  ^(A',^')  satisfying 


a) 

* "■'d-j+1 

(10.8) 

b) 

= -1  and 

c) 
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The  requirements  (10.8)(a)  and  (10.8)(c)  are  compatible  since  the 

number  d of  basis  functions  is  even.  Let  9 - F.,v  * F , , v be  the 

A A ,m 

quasilnterpolant  of  v with  points  {Xj}.  At  this  point,  all  we  know  is 
that  F^.v  « s'‘(A'  ,z'). 


By  the  symmetry  of  the  function  v,  the  partition  and  the 
quasiinterpolation  points  ^ ^ » 


(v)  = X (v)  , 1 < j < d. 

J d-j+1  - - 


Thus , 


Z A.(v)  N.  , 

j = l J J’*" 

d/2 

jf,  <“j,k  * Vj+i.k> 

d/2 

Z X.  B. , 


j = l 


J J 


which  shows  that  F.,v  ^ ES  (A',z').  By  (10.8)(b)  and  Lemma  3.3,  F -v 
A — A 

interpolates  v at  -1  and  1.  Thus,  F^,v  e ESq(A',_z')  and 


9 = F^,v  c SS^(A,z). 

[0,1] 


We  now  consider  error  bounds.  Since  we  are  concerned  with 
restrictions  to  [0,1],  we  need  to  bound 

X . 


x”  ^(D*^(v  - F^v))^  dx,  1 ^ i < N. 


1-1 


We  consider  two  separate  cases,  depending  on  whether  or  not 

®1  ® (-X^,Xj^)  ■»  (?. 


12A 


Let  1 be  such  that  9^  9 (-XpX^  4 (?.  Then  0 < t £ k.  Let  b be 
the  smallest  number  such  that  [~b,b)  and  b ^ x^*  Then 


(10.9) 

Moreover, 

(10. 10) 


h,  < b < X.,,  , < (2k-l)h. 
1 — i+k-1  — 


For  all  1 £ j £ 1> 
2k-l 

^ ^ if! 

min  h, 
l<l<i  ' 


2k-l 

£ I (M(A))*' 
1-1 

< (M(M)^^-l  = 

- M(A)-1 


(see  Figures  10.1  and  3.2). 


Let  S V be  the  local  polynomial  approximation  to  v given  by  Lemma 
b 

10.2,  and  extend  S v to  (-b,b)  by  reflection  about  the  y-axis: 
b 


Sj^(-x)  - Sj^(x),  0 £ X £ b. 
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i)  is  a polynomial  of  degree  m-1  in  each  of  the 

intervals  (-b,0)  and  (0,b); 


11)  is  even; 


iii)  S,  V has  in  + 1 continuous  derivatives  at  0. 

D 


For  (iii),  note  that  the  (one-sided)  even  derivatives  agree  at  0 since 

S V is  symmetric  about  the  y-axis,  while  the  odd  derivatives  of 
b 

appropriate  orders  all  vanish  at  0 by  (10.4). 


Let  S V be  defined  for  x i (-b,b)  by  extending  each  of  its  two 
b 

Ir 

polynomial  pieces.  Clearly  S^^v  e ES  ,z').  Let  E^^v  = v - S^^v.  Since 
the  quasiinterpolant  reproduces  elements  of  ES 


L»(l^) 


(10. 11) 


We  want  to  get  bounds  like  those  of  (10.3)  for  E,  . In  view  of 

l.t 

(10.11),  we  only  need  to  bound  II  , 2 n*  By  (3*^) 

A b L (Ij^) 

and  the  triangle  Inequality, 
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By  Lenuna  3.4,  (10.8)(c),  and  (10.9), 


|A^(Ej^v)|  < Z |d‘'e^v(Tj)| 


< C(k,M(A))  hj 


< C(k,M(A))  ||D*^E^v 


b " L”(x^,b) 


1 1 ^4 

Since  b > x„,  — r < — < Thus,  by  Lemma  9.4,  Lemma  10.2,  and 

- 2 (b  - X,  ) - - b ’ ^ 


r 2 


(10. 10), 


lXj(Ej^v)|  < C(k,M(A))  b"^  Cj  x[""  [(b-xp“^  llo’^E^v  II 


r ( \r  r*  « « (l*n)/2  t m“  r*  ( 1 /2  ) ■ ■ _^m  • ■ 

£ C(k,M(A))  b Cj^Cj^c^  x^  b ||  D v || 

(When  Tj  < 0,  we  use  the  fact  that  E^v  is  even,  whence  (D*^Egv(T^)(  = 
iD’^EgV(-Tj)  I .)  Thus,  by  (10.  10), 

|X^(EgV)|  < C(k,M(A))C^c^  ll°”''llB(b) 

Together  with  (3.2)(b),  Lemma  10.3,  and  (10.10), 


127 


II  (n-1  )/2  ,1 

II  X D 


< ^ |\(E.v)|  llx^""^^''^  D*'n 


jfv. 


j'  b 


j.k  "l*(I^) 


(10. 12) 


B(b) 


< kc2C3c;b”-^  ll'^‘“v||3(,j. 


It  follows  from  (10.11),  (10.3),  (10.12),  and  (10.9)  that 


^ ^*=1  kc2C3cJ]  b““^  llD”vllB(b) 


(10. 13) 


< [c^  + k.C2C3C.Jl  (2(k-l))'^^  h®"*'  ||D®v|| 


B(b) 


ll°"''”B(b)- 


We  now  turn  to  the  case  that  9.  n (-XpX^  = Ql.  Let  x,  be  the 

^ i 

first  knot  to  the  left  of  9^  (see  Figure  6.1).  Note  that 


I 6^ I = length  of  9^  < 


2kh 


and 


— - < k. 
■^i 


We  have,  using  the  error  bounds  of  Theorem  3.1, 


[ 

i 

! 
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*1 

* X ^ 

i-1 

< (K(2k)"*"‘)2  h2(®-‘)  / (D“v)2 

®i 

< (K(2k)“-S2  ( , jn-1  ^2(m-l)  ^ 

^ Ji  e^ 

< (K(2k)”"*’)^  k""^  j^2(iii-4)  j.  ^n-l^jjm^j2 


c2  / x"-^d“v)2  dx. 


Together  with  (10.13)  and  the  usual  observation  that  not  more  that 
2k  of  the  0^  intersect  at  any  point,  we  have 


d\v  - F,.v)  II  2 


2 2 2 
where  = 2k  max(c^,  c^). 


y 2 ^ ^2  2(m-J,)  ,,  m . 2 


D 


We  now  consider  the  error  in  the  RRG  approximation  to  the 
generalized  solution  u of  (9.4)  - (9.5)  and  show  that  the  RRG 

mm  k 

approximation  u in  SSq(a,^)  is  an  optimal-order  approximation  to  u. 
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Theorem  10»2;  Let  u be  the  generalized  solution  of  (9.4)  - (9.5).  let 
u c SSq(a,£)  be  the  RRG  approximation  to  u.  If  f e then 


Hu  - u llg  < (AC3r)^  h^  Ilf  II 

where  i = min(m,  k) . 


Proof:  By  Lemma  9.2,  the  generalized  solution  u g J®  n The  error 

bounds  then  follow  from  the  approximation  results  of  Theorem  10.1,  the 
regularity  theorem  (Lemma  9.2),  and  Nitsche's  trick  (Theorem  2.6). 


□ 


10. 4 Computational  Aspects  of  the  Method 

Other  authors  who  have  treated  the  problem  (9.4)  - (9.5)  have 
imposed  either  no  constraint,  or,  at  most,  the  natural  boundary 
condition  Du(0)  = 0 on  the  space  of  approximate  solutions  ( [R3] , [J3], 
[D6]).  The  resulting  approximation  to  u,  while  accurate,  will  not  give 
rise  to  a smooth  function  in  the  original  domain  B c R*'.  By  using  the 
space  SSq(a,£),  this  smoothness  is  obtained. 

An  alternative  is  to  deal  with  the  problem  (9.4')  - (9.5').  In 
that  case,  rather  than  enforcing  the  requirement  of  evenness  on  the 
subspace,  as  we  did  with  ESq(a',£')  (and  SSq(a,£)),  we  take  the  RRG 
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v 

approximation  to  u from  the  space  ,z')  . Since  the  basis  functions 
are  no  longer  even,  we  would  have  to  take  the  Inner  products  over  (-1,1) 
instead  of  (0,1).  Surprisingly,  the  restriction  of  this  spline  to  [0,1] 
is  the  RRG  approximation  to  u from  SS^(A 


It  suffices  to  show  that  the  RRG  approximation  u c 

even.  Clearly,  by  the  symmetry  property  (10.1)  of  the  B-splines, 
d 

u ^ ,k  even  if  the  vector  £ is  symmetric  about  its  middle, 

i.e. , 


(10.14)  Cj  - 1 < j < d. 

These  coefficients  are  obtained  as  the  solution  of  the  linear  system 
Ax  » ^ of  (2.10).  Because  of  the  evenness  of  the  data  (functions  p,  q, 
and  f)  and  the  basis  functions,  the  matrix  A will  be  symmetric  about 
the  alternate  diagonal  and  the  vector  f^  will  be  symmetric  about  its 
middle.  Thus,  the  coefficients  of  u will  satisfy  (10.14). 

It  might  appear  that  by  dealing  with  a two-point  boundary  value 
problem  on  (-1,1)  in  which  all  the  functions  involved  are  even,  we  would 
do  twice  as  much  work  as  is  necessary.  This  is  Incorrect.  It  does  not 
cost  any  more,  in  work  and  storage,  to  use  the  (-1,1)  problem  than  the 

If 

(0,1)  problem.  Suppose,  for  example,  that  we  use  the  spac  ‘^j(d',£')  . 
Because  of  the  symmetries  of  Che  matrix  A,  only  1/4  of  its  elements  need 
CO  be  computed.  Moreover,  using  an  algorithm  of  Evans  and  Hatzopoulos 
[E3]  which  Cakes  advantage  of  symmetry  about  Che  alternate  diagonal,  Che 
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equations  (2.10)  can  be  solved  In  half  Che  time  required  by  Che  usual 
band  Cholesky  algorithm. 

The  effect  of  numerical  quadrature  (used  to  compute  the  matrix  A 
and  the  right-hand  side  vector  f^)  on  the  accuracy  of  the  RRG 
approximation  has  been  analyzed  by  Fix  for  nonsingular  problems  [FI]. 

He  showed  that  If  the  integrals  are  computed  using  composite  Gaussian 
quadrature  with  k-1  points  in  each  interval,  then  the  error  due  to  the 
quadrature  is  asymptotically  as  small  as  the  discretization  error.  We 
conjecture  that  this  result  applies  to  the  singular  problem 
(9.4)  - (9.5)  and  that  k-1  points  suffice.  The  numerical  results  of  the 
next  section  strongly  support  this  viewpoint. 


10.5  Numerical  Resul ts 

In  this  section  we  present  the  results  of  a numerical  experiment, 
which  illustrates  the  utility  of  the  computational  procedure  analyzed  in 
the  previous  sections.  Following  Russell  and  Shampine  [R3] , we  consider 
the  problem 

- D(x^Du)  + 4x^u  = - 20x^ 

u(-l)  = u(l)  = 0 

which  has  the  solution  u(x)  » ^_sinh^  _ 

X sinh  2 

The  RRG  approximations  to  u from  several  of  the  ES^(4'^2')  spaces 
were  computed,  and  the  error  tabulated  below.  The  partition  4'  of 
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(-1,1)  was  uniform,  with  2N  subintervals  and  mesh-width  h - All 

N 

computations  were  performed  in  double  precision  on  a PDP-10  (with  5A 
binary  digits).  The  integrals  required  were  computed  using  Gaussian 
quadrature,  with  k-1  nodes  in  each  interval  of  the  mesh.  We  give  the 
norms  II  e H ^ and  ||  De  ||  g,  computed  using  k+l  - point  composite  Gaussian 
quadrature  rules,  and  also  the  quantity  II  Further  details 

concerning  notation  and  the  observed  rate  of  convergence  (RATE,  below) 
are  given  in  Section  6.4. 

As  predicted  by  the  theory,  the  rate  of  convergence  appears  to  be 
Ic  VI 

h for  the  error  and  h*^"^  for  the  derivative.  Moreover,  k-1  quadrature 
nodes  per  interval  are  sufficient  to  maintain  the  predicted  rate  of 
convergence. 


(-01) 

12  (-01) 

1.98 

[ 

54  (-02) 

2.00 

31  (-02) 

2.00 

1 

20  (-02) 

1 2.00 

14  (-02) 

2.00 

10  (-02) 

2.00 

76  (-03) 

2.00 

N A 


.38  (+00) 
.12  (+00) 
.58  (-01) 
,35  (-01) 
.23  (-01) 
,17  (-01) 
,13  (-01) 
.10  (-01) 


.42  (+00) 
,21  (+00) 
,14  (+00) 
.10  (+00) 
.82  (-01) 
.69  (-01) 
,59  (-01) 
,51  (-01) 


Table  10.1:  Error  in  ESq(a') 


rate 


,11  (-02) 
.10  (-03) 
,28  (-04) 
.11  (-04) 
.57  (-05) 
.33  (-05) 
.21  (-05) 
.14  (-05) 


15 

(-02) 

12 

(-03) 

26 

(-04) 

90 

(-05) 

40 

(-05) 

21 

(-05) 

12 

(-05) 

73 

(-06) 

.19  (-01) 

.46  (-02) 

2. 

.20  (-02) 

2. 

.11  (-02) 

2. 

.72  (-03) 

2. 

.50  (-03) 

2. 

.37  (-03) 

2. 

.28  (-03) 

2. 

Table  10.2:  Error  in  ES^(A') 


r ' ^ 
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N 

He 

N II  B 

RATE 

He 

N II  «,  A 

RATE 

II  II  B 

RATE 

■ 

.41 

(-04) 

.96 

(-04) 

.11 

(-02) 

B 

.28 

(-05) 

3.83 

.64 

(-05) 

3.91 

.15 

(-03) 

2.86 

12 

.59 

(-06) 

3.86 

.14 

(-05) 

3.72 

.46 

(-04) 

2.90 

16 

.19 

(-06) 

3.90 

.47 

(-06) 

3.80 

.20 

(-04) 

2.93 

20 

.81 

(-07) 

3.93 

.20 

(-06) 

3.84 

.10 

(-04) 

2.94 

24 

.39 

(-07) 

3.94 

.99 

(-07) 

3.87 

.60 

(-05) 

2.95 

28 

.21 

(-07) 

3.95 

.54 

(-07) 

3.89 

.38 

(-05) 

2.96 

32 

.13 

(-07) 

3.96 

.32 

(-07) 

3.90 

.26 

(-05) 

2.96 

Table 

10.3: 

Error 

in  ESq(A 

') 

N 

He 

N II  B 

RATE 

He 

N OD,  i 

RATE 

11  H b 

RATE 

m 

.18 

(-06) 

. 12 

(-05) 

.42 

(-05) 

.36 

(-08) 

5.65 

.82 

(-08) 

7.21 

.17 

(-06) 

4.60 

.33 

(-09) 

5.91 

.76 

(-09) 

5.85 

.24 

(-07) 

4.85 

.60 

(-10) 

5.95 

.14 

(-09) 

5.86 

.59 

(-08) 

4.91 

20 

.16 

(-10) 

5.97 

.38 

(-10) 

5.88 

.19 

(-08) 

4.94 

24 

.53 

(-11) 

5.97 

.13 

(-10) 

5.89 

.79 

(-09) 

4.95 

28 

.21 

(-11) 

5.98 

.52 

(-11) 

5.90 

.37 

(-09) 

4.96 

32 

.96 

(-12) 

5.97 



.24 

(-11) 

5.89 

. 19 

(-09) 

4.97 

Table  10.4:  Error  in  ES^(A') 


CONCLUSIONS 


i 
) 

j 

i 

3 

I 

For  the  two  important  classes  of  singular  two-point  boundary  value 
problem,  we  have  shown  how  finite  element  methods  can  accurately 
approximate  the  solution.  We  feel  that  for  linear,  one-dimensional 
problems  with  a singularity  at  one  of  the  endpoints,  many  of  the 
important  problems  have  been  solved.  Also,  our  results  probably  can  be 
extended  without  great  difficulty  to  time-dependent  and  mildly  nonlinear 
problems. 

I In  a few  areas,  we  have  no  satisfactory  theoretical  results.  How 

I much  B-gradlng  is  required  for  the  subspaces  of  Chapter  6?  It  seems 

2 

j that  a 62-graded  mesh  is  the  thing  to  use  (if  L -norm  error  is 

I 

I important) , but  we  have  not  been  able  to  prove  this.  The  lack  of  any 

' local  error  bound  for  the  RRG  approximation  has  defeated  our  attempts. 

I • 

We  have  said  nothing  about  the  collocation  method.  Nevertheless, 

theoretical  results  for  nonsingular  problems  and  numerical  experiments 

for  singular  problems  indicate  that  it  is  an  attractive  alternative  to 

j RRG  for  singular  two-point  boundary  value  problems.  Using  splines  on  a 

I B-graded  mesh  (as  in  Chapter  6),  collocation  produces  an  optimally 

‘ 2 

I accurate  L -norm  approximation,  provided  we  collocate  at  the  Gaussian 

t 

points.  Unfortunately,  these  methods  apply  only  to  C^  spline  spaces. 

(, 

L 
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The  results  of  a numerical  experiment  on  the  problem  of  Section  6.4 

are  given  In  Tables  C. 1 - C.3.  It  appears  that,  unlike  RRG,  collocation 

2 

succeeds  In  producing  an  optimal-order  L -norm  approximation  In  the 

spaces  ^)i  Other  experiments  show  that  collocation  at  the 

weighted  Gaussian  points  Is  effective  In  the  weighted  spline  spaces  of 
Chapter  7. 


We  have  not  given  any  theoretical  justification  for  the  use  of 
numerical  quadrature  in  the  RRG  method.  It  appears  that,  just  as  in  the 
nonsingular  case,  k-1  quadrature  nodes  In  each  interval  are  sufficient 
to  preserve  the  rate  of  convergence,  but  we  have  no  proof  of  this. 


The  generalized  L-spline  spaces  of  Chapter  8 suffer  from  a severe 
restriction:  beyond  continuity,  no  smoothness  my  be  imposed.  For  the 
same  partition,  a smooth  space  would  have  smaller  dimension,  making  the 
RRG  approximation  less  costly  to  compute. 


For  the  spherically  symmetric  problem,  we  feel  that  the  theory  is 
more  solid.  Jesperson  has  already  succeeded  in  obtaining  L°°-error 
bounds  IJ3],  while  de  Hoog  and  Weiss  have  some  results  on  collocation 
[D4].  The  effect  of  numerical  quadrature,  however,  still  remains  to  be 
analyzed. 


Table  C.2:  Collocation:  Error  in  S 


8 

.36  (-04) 

.12  (-02) 

.40  (-01) 

16 

.69  (-06) 

5.69 

.15  (-03) 

3.00 

.14  (-01) 

1.50 

24 

.69  (-07) 

5.69 

.45  (-04) 

3.00 

.77  (-02) 

1.50 

32 

.13  (-07) 

5.72 

.19  (-04) 

3.00 

.50  (-02) 

1.50 

48 

.13  (-08) 

5.79 

.56  (-05) 

3.00 

.27  (-02) 

1.50 

Table  C.3:  Collocation:  Error  in  5 e 
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